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Abstract 

Topological invariants built from the periodic Bloch functions characterize new phases of matter, 
such as topological insulators and topological superconductors. The most important topological 
invariant is the Chern number that explains the quantized conductance of the quantum Hall effect. 
Here, we provide a general result for the superfluid weight Dg of a multiband superconductor that 
is applicable to topologically nontrivial bands with nonzero Chern number C. We hnd that the 
integral over the Brillouin zone of the quantum metric, an invariant calculated from the Bloch 
functions, gives the superfluid weight in a flat band, with the bound Dg > \C\. Thus, even a flat 
band can carry hnite superfluid current, provided the Chern number is nonzero. As an example, 
we provide Dg for the time-reversal invariant attractive Harper-Hubbard model that can be exper¬ 
imentally tested in ultracold gases. In general, our results establish that a topologically nontrivial 
flat band is a promising concept for increasing the critical temperature of the superconducting 
transition. 


I. INTRODUCTION 

An important result of Bardeen-Cooper-Schrieffer (BCS) theory is the relation oc 
exp (^~ uno ’ between the critical temperature of the superconducting transition and 
the microscopic parameters of a superconductor, such as the coupling constant U of the 
effective attractive interaction and the density of states at the Fermi energy no (hip)- This 
result is valid in the limit where the coupling constant U is much smaller than the band¬ 
width, which is roughly given in a tight-binding approximation by the hopping energy J 
between neighbouring atomic orbitals. The BCS formula suggests two ways to increase the 
critical temperature, namely either to enhance the coupling constant U or the density of 
states no (Tip). Whereas the electron-electron attraction parametrized by U is the result 
of complicated many-body physics, not yet well-understood in the case of unconventional 
superconductors, the density of states can be more easily obtained and engineered in a 
single-particle framework by means of band structure calculations. 

The density of states at the Fermi energy no(T^p) is maximal for vanishing bandwidth 
and so is the critical temperature. In this limit, the energy dispersion as a function of lattice 
quasi-momentum Kk is constant e(k) = d and the corresponding energy band is called a “flat 
band”. The exponential suppression of the critical temperature disappears in the flat-band 
limit f//J ^ 1 since BCS theory predicts^^^ oc Uno{Ep) oc U/J. This might provide the 
way to reach the grand goal of room-temperature superconductivity. 

A crucial question unaddressed in many works on flat-band superconductors^^^’^^^^ is 
whether the superfluid mass density ps, or better, superfluid weight Dg (see below), is 
nonzero, leading to the Meissner effect and dissipationless transport^^’^^ that dehne su¬ 
perconductivity. Within the single-band effective Hamiltonian approximation^^^^®, in which 
only the band dispersion enters, the superfluid weight vanishes {Dg oc J) since Cooper pairs 
localize in the individual lattice sites. Finite superfluid currents can be found in some hat¬ 
band systems^^, but a general theory, connecting the superhuid weight to invariants of the 
band structure (possibly topological invariants) has not yet been provided. The aim of this 
work is to answer, at a general level, the crucial question whether superhuidity can exist in 
a hat band and to explore its possible connections with topological properties of the band. 

Using a multiband BCS framework, we show that the superhuid density depends not 
only on the energy dispersion but also on the Bloch functions of a lattice Hamiltonian. This 
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fact is especially important in the flat-band limit. Moreover, we argue that the superfluid 
density is subtly affected by the topological invariants encoded in the Bloch functions even in 
conventional superconductors (not topological)^®. Topological invariants such as the Chern 
number C are gauge-invariant integer-valued quantities^®’^^ which determine the charge and 
spin conductance and the presence of robust edge states^^^^®. Indeed, the physical picture 
of localized Cooper pairs is intimately related to the existence of exponentially localized 
Wannier functions^^ that can be constructed only if the Chern number C is nonzero^® (see 
Fig. 1). Note that the Chern number corresponds to an antisymmetric tensor, the Hall 
conductance, whereas the superfluid weight is a symmetric one and, if nonzero in a flat 
band, is an invariant quantity constructed only from the Bloch functions. We hnd that 
the superfluid density in a flat band is proportional to a symmetric tensor given by the 
Brillouin-zone average of a quantity known as the quantum metric^®’^^. This tensor is the 
real part of an invariant matrix M., which depends only on the Bloch functions, while the 
imaginary (antisymmetric) part is the Chern number. By means of the properties of the 
invariant M., we prove a bound on the superfluid weight that reads Dg ^ |C| in appropriate 
units (see Fig. 1). Moreover, we predict that the superfluid weight is proportional to the 
coupling constant Dg oc t/ in a flat band. As a concrete application, we derive the superfluid 
weight in closed form for the Harper-Hubbard modeF®. Using artihcial gauge helds, the 
Harper model has been recently realized with ultracold gases^®’®°, which are a good platform 
to verify our predictions. Our arguments are general and similar results are expected for 
other flat bands or bands that are only partially flat. 


II. RESULTS 


A. Effective lattice Hamiltonian 


Our goal is to provide, within a mean-held approximation, a general formula for the 
superhuid weight of a multiband system which can include topologically nontrivial bands 
and/or hat bands. A hnite supercurrent is associated with a winding of the phase of the 
superconductor complex order parameter A(r). In the specihc case of a constant current 
J(q). the order parameter has the form of a plane wave A(r) = |A|e^*'^''' with wavevector 2q. 
The superhuid mass density ps and superhuid weight are dehned as the change in the free 
energy density AF/V = ^PsV^ = |T*sPs (U is the volume in three dimensions, or the area in 
two dimensions), due to the motion of Cooper pairs with uniform velocity Ug = h\q\/m and 
momentum ps = 2h|q|. In lattice systems, the mass m is not a well-dehned concept, and it 
is better to use the superhuid weight^^’^® Dg. 

A computationally convenient dehnition of the superhuid weight is in terms of the grand 
potential f2(T, p, A,q) (see Ref. 31 and Supplementary Note 1) 


[B.i 




1 


Vh? dqidqj 


/L,A,q=0 


( 1 ) 


where i, j = x,y,z are spatial indices. In anisotropic and time-reversal invariant systems, the 
superhuid weight is given by a symmetric tensor [Ds\ij (the notation [M]ij for the elements 
of a matrix M, with i,j not necessarily spatial indices, is used throughout the article). 

In calculating the superhuid weight, we proceed in the following way. 1) The supercurrent 
wavevector q is introduced in the Hamiltonian in a way that is rigorous for topologically 
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non-trivial bands: a multiband approach is used. 2) The kinetic Hamiltonian is Fourier 
transformed, which dehnes the band dispersions and the Bloch functions. 3) A mean-field 
approximation is done by introducing a Bogoliubov-de Gennes (BdG) Hamiltonian. 4) The 
BdG Hamiltonian is diagonalized to provide a convenient expression for the grand potential. 

5) Supercurrent and superfluid weight are obtained as derivatives of the grand potential with 
respect to q: the results are given in terms of the band dispersions and the Bloch functions. 

6) The results are connected to topological properties of the system. 

Some care is needed to introduce the wavevector q in the Hamiltonian in a proper way. 
By a suitable gauge transformation, it is possible to constrain the complex order parameter 
A(r) to be real and have the same translational symmetry as the underlying lattice, whereas 
the wavevector q appears in the kinetic term of the lattice Hamiltonian 

^ , ( 2 ) 

where the matrix elements iFi j oc J are hopping amplitudes between lattice sites. If the 
wavevector q is identified with a constant external vector potential A according to q = gA/h, 
Eq. (2) becomes the usual Peierls substitution. 

The Peierls substitution is an approximation valid only if the basis states of the lat¬ 
tice Hamiltonian are well-localized^^’^®, and ideally they should be exponentially localized 
Wannier functions^^. Since bands with a nonzero Ghern number do not allow exponentially 
localized Wannier functions^®, we use a multiband approach that can circumvent this prob¬ 
lem. We consider a subset of bands, which we call S, well-separated from other bands by 
band gaps (a composite band)^^ such that the Ghern number (or numbers) of the composite 
band is zero. By linear superposition of Bloch functions of all the bands in S, it is possible 
to construct exponentially localized Wannier functions. For the notation and the definition 
of Wannier functions, see Figure 2 and Supplementary Note 2. 

In the basis of Wannier functions, the effective lattice Hamiltonian for the composite 
band reads (the derivation is known^^, but for convenience we repeat it in Supplementary 
Note 2) 


it- = 

iaj/3 a 




^ - /iA . 


( 3 ) 


The are annihilation (creation) operators for the orbitals (Wannier functions) la¬ 

belled by \a and j/3 (Fig. 2a) and spin a, fi the chemical potential, N = X^iao-the 
particle number operator, and we consider the specihc case of an attractive Hubbard interac¬ 
tion [U >0). The Peierls substitution has been used, properly generalized to the multiband 
case (see Fig. 2a for the definition of ria). For q = 0, the Hamiltonian is invariant under 
time-reversal symmetry (TRS) since (Aj|^j^)* = and invariant under spin rotation 

around the 2 :-axis, but in general ^ K^. 

Diagonalization of the Fourier transform of the hopping matrix in Eq. (3) gives the band 
structure (see Methods IV A) 

(4) 

Here = diag(e„k( 7 ) is a diagonal matrix composed of the dispersions Snka of each band (n 
labels a single band belonging to the composite band), while the n-th column of the unitary 
matrix ^ko- is the Bloch function [^ko-]a,n = fi'nko-(tt) of the n-th band. TRS implies that 
£kt = ^-k; and = G-ki- 
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B. BCS theory and superfluid weight in a multiband system 


The idea of superconductivity in multiband systems dates back to^^ 1959. The hrst 
superconductor for which multiband effects are indeed measurable is magnesium diboride 
(MgB 2 ), discovered as recently as^^’^^ 2001. However to the best of our knowledge, a general 
and consistent theory for the superfluid weight in a multiband system, in particular for 
topologically nontrivial flat bands, has not yet been worked out. 

In the following, we develop the theory of the superfluid weight in a multiband system 
within the framework of BCS theory, namely, we use a mean-held decoupling of the inter¬ 
action term 




+ H-C. 


( 5 ) 


where AiQ, = —U{ciaiCia^). Furthermore, we choose a gauge where the gap function preserves 
the discrete translational symmetry of the lattice, which means Ai^ = A^. Thus the pairing 
terms are diagonal in momentum space Xlia ~ Xlka already 

discussed, the wavevector q enters in the mean-held Hamiltonian in the kinetic energy term 
through the Peierls substitution (3). 

In terms of new fermionic operators d^ko- (see Methods IV A) the mean-held Hamiltonian 

reads 'Hm.f. = Xlk ^k-^k(q)dk with the Nambu spinor given by dk = (dnkt) the 

BdG Hamiltonian is a 2 x 2 block matrix dehned by 


^k-q hi ^k_qA^k+q ^ 

^k+q^^k-q ~ (^k+q ~ hi) j 


( 6 ) 


with A = A* = diag(AQ,). Due to TRS, the spin index a has been dropped, and only 
the eigenvalues and eigenstates for the spin up are used in the following, namely £k = ^kt 
and ^k = Q\c\- The BdG Hamiltonian is diagonalized in terms of a diagonal matrix of 
quasiparticle excitation energies F'k(q) and a unitary matrix VVk(q) 

HM) = Wk(q)i?k(q)W^(q) . (7) 


The symmetries of the BdG Hamiltonian for q = 0 imply that these matrices have the 
following structure (i^nk > 0, see Supplementary Note 3) 


B.(q = 0) = 


Wk(q = 0) = 


Wk 

Vk 


0 

-diag(.F„k; 

-Vk 
Wk 


( 8 ) 

(9) 


While the kinetic energy terms of the BdG Hamiltonian (6) are diagonal in the band index, 
the pairing terms depend in a complicated way on the Bloch functions and on the order 
parameters A^ relative to all orbitals. It is interesting to explore the consequences of this 
nontrivial structure on superfluid transport. A “gauge” transformation of the Bloch func¬ 
tions given by ^k —t ^kx4k, with Mk a unitary matrix subject to the constraint of commuting 
with the matrix of band dispersions [^k, x4k] = 0, leaves the BdG Hamiltonian (6) unchanged 
in form while the eigenfunctions (9) change accordingly. This freedom in the dehnition of the 
Bloch function is the same one preventing a unique dehnition of Wannier functions^^. All 
observable quantities, such as current and superfluid weight, are necessarily gauge invariant. 
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At zero temperature, the grand potential is (see Methods IV B and Supplementary Note 

3) 

fi(q) = -^ETt|ffk(q)ll + .... (10) 

k 

The dots in the above equation represent terms in the grand potential that do not contribute 
to the superfluid weight. The superfluid current density is obtained from the first derivative 
of fl(q) 


J(q) 


-1 

Wh 


^Tr sign(T;k(q))>Vj(q)(9qi/k(q)>Vk(q) 


with — (9qhfk(q) 


( -^k^k-q <9qT>k(q)\ 

V-'9q^^k(-q) l^^k+q ) 


( 11 ) 


The definition T’k(q) = — ^k_qA^k+q has been employed above. Due to the linearity of 
the trace in Eq. (11), the current splits into two contributions that are separately gauge 
invariant. We call the first the “conventional” current, which depends on the group velocity 
(9k£k/h and is of order J/h, and another contribution of order A/h that comes from the 
off-diagonal blocks in Eq. (11). Our prediction of the latter current component is highly 
interesting since it may be nonzero in a flat band, unlike the conventional component. Note 
that in the semiclassical expression for the velocity in a magnetic Bloch band^^’^®, two terms 
appear as well; the group velocity obtained by the band dispersion and the Berry curvature, 
which is due to interband coupling as the off-diagonal blocks in Eq. (11). However, the 
analogy is not complete and the precise relation between Berry curvature and the interband 
contribution to the superfluid density is clarified below. 

The superfluid weight is obtained by taking the derivative of the current density J(q) 
and setting q = 0. The superfluid weight consists of three terms Dg = i -|- Ds ^2 + 
(details of the derivation are provided in Supplementary Note 3). We call the first term the 
conventional superfluid weight 




vh? ^ 

k 


Tr 




( 12 ) 


This is the only term present in the single band case, and is zero for a flat band. The other 
terms are present only in the multiband case. The second term stems from the derivative 
^q^ of the off-diagonal blocks in Eq. (11) 




Vh^ ^ 

k 


Tr 


VkWi^a,,9,^Dk(q = 0) 


(13) 


Finally, we have a contribution from terms of the form VVj^(q)(9qVVk(q); 


[Ds,3\i,j - y^2 X] X] 


ri,n' [-^k, 




En'k + E, 


n'k 


(14) 


Bk, = UldqPM = 0)Wk + VidqPM = 0)Vk + Vldu^e^U^ - ul^k,e^v ^,. 

Eqs. (12)-(13)-(14) are the main result of our work since the superfluid weight can be readily 
calculated using only the ground state solution (8)-(9). The conventional superfluid weight 
Dg 1 is invariant under gauge transformations, which means that Dg 2 + Dg 3 is itself gauge 
invariant, thus the superfluid weight splits into two distinct contributions in the same way 
as the current. 
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C. Superfluid weight in a flat band 


The general results in Eq. (12)-(14) can be specialized to the case of a flat band in two 
dimensions, and a particularly interesting case is that of a topologically nontrivial flat band. 
A band specihed by n within the composite band S is considered for which the band gaps 
separating it from the lower (h — 1) and upper (h -|- 1) bands are large with respect to the 
bandwidth. It is thus possible to have a coupling constant U such that 


mink£(n+i)k - maXkSfik 
miuk^nk - maxk£(n-i)k 


3> 3> maxk^fik — mink^nk • 


(15) 


In this limit, the dispersion of the h-th band can be approximated by its average Enk ~ 

To proceed, it is assumed that the order parameters Aq, relative to each orbital in the unit 
cell are all equal, in other words, that T>k(q = 0) = Al where A = A^ is now a real scalar. 
We can prove this fact rigorously for the Harper-Hubbard model. 

Given this assumption, it is shown in Supplementary Note 4 that an approximate self- 
consistent solution can be found in the limit (15) and has the following form. The matrices 
[Wk]n,n' = Un5n,n' and [Vk]n,n' = VnSn,n' are diagonal, while the other relevant quantities are 


Ur,. = 


i + 


(16) 

mUfi U Tlfj) 

\/z/(l - v) , 

(17) 

0 < n < h 



n = n 

Vn = 

(18) 

n > n 



jp U Tl(j) 

— 2 ’ 

I'^nk /^| 7 

n = n 

n ^ n, 

(19) 

and = . 

Norh (the number of orbitals; 

see Fig. 2a). 


Ert.U — 


This solution depends only on Eq. (15) and is in fact generic for any flat band. The only 
assumption is A^ = A. 

The above solution can be inserted in the general formulas (12)-(13)-(14). The conven¬ 
tional superfluid weight Dg^i vanishes in the flat band limit, and the remaining part has the 
form 

2Und 

- 7/1 \ — 7 /) A/,.. 


= [Ds,2h + . ( 20 ) 


Tlh? 


We thus And in the flat-band limit that the superfluid weight is proportional to A oc Un^. 
This is consistent with Ref. 17 for the specific case of the flat band of surface states in rhom- 
bohedral graphite, however, our theory is much more general and can be applied to a variety 
of systems. This result has to be contrasted with the one for an ordinary superconductor 
in a parabolic band Dg = n^/meEi oc J (with n-p the total particle density and rriefr the 
effective mass) that can be obtained from Eq. (12), the only term that survives in the single 
band case. Therefore, an important prediction is that in a flat-band superconductor, the 
superffuid weight is linearly dependent on the coupling constant, whereas it is independent 
from it in an ordinary superconductor. Interestingly also in superconducting graphene with 
the chemical potential tuned at the Dirac point, one has^^’^® Dg oc U. 
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(21) 


The matrix = Re(Aljj) is the real part of a Hermitian matrix dehned as 

= d^kB,,{k) 

Jb.z. 

which is the integral over the whole Brillouin zone of the so-called quantum geometric tensor 

B,,{k) = 2Tr [{dk,gl){dk^g^)'\ + 2Tr , (22) 

where is the projection of ^k on the h-th band (see Methods IV C). In mathematical 
terms, the quantum geometric tensor is the Fubini-Study metric in the projective manifold of 
quantum states^®’^’^. The quantum geometric tensor has been recently related to observable 
quantities in a single-particle context such as the noise current spectrum^®, and plays an 
important role in characterizing bands that can host fractional Chern insulators, namely, 
lattice generalization of the fractional quantum Hall state^°’^^. 

It can be shown that Bij (k) is zero if ^k is a square unitary matrix. The case where ^k is 
a square matrix corresponds to superfluid pairing including all the bands of the composite 
band. Consistently, only if a strict subset of the bands in the composite band participate in 
the pairing, the superfluid weight can be nonzero in the flat-band limit. In contrast in the 
same limit, the whole composite band, which has zero Chern number, is a set of localized 
orbitals with vanishing hopping. The imaginary part of Bij{k) is the well-known Berry 
curvature, and its integral over the Brillouin zone is the Chern number in two dimensions 
Im(Aljj) = Mjj = €ijC (cij = —€ji is the Levi-Civita symbol). The Chern number refers to 
spin-resolved bands since the component of the spin is a conserved quantity. 

The real part of Bij(k) is a Riemannian metric^®’^^ dehned over the Brillouin zone, the 
so-called quantum metric. In two dimensions, the positive semidehniteness of the 2x2 
matrix Mij (see Methods IV C and Fig. 1) implies that det(Al) > 0 or 

det(Af^) > det(Af^) = . (23) 

For an isotropic system, the matrix is proportional to the identity, and this gives the 
bound Dg > |C|, in appropriate units. The trace TrAl is the gauge invariant part of 
the localization functional for Wannier functions F studied by Marzari and Vanderbilt^®’^^, 
pointing to an intimate connection between non-localization of Wannier functions and the 
superhuid weight. Indeed, Eq. (23) also implies that the localization functional is bounded 
from below by the Chern number. In 2D, the bound is F > with the area of the 

unit cell (see Supplementary Note 5). 


D. Time-reversal invariant attractive Harper-Hubbard model 


To make our results more concrete and study the superhuid weight in a quasi-hat band, 
we consider the specihc example of the time-reversal invariant attractive Harper-Hubbard 
modeF®. This model is dehned on a two dimensional square lattice with lattice spacing a 
by the hopping operator 


Kfj = -J(u,- 




with X = (1, 0)"^, y = (0, l)'^ and lo = The phase factors are the lattice version 

of the Landau gauge that introduces a uniform magnetic held with hux per plaquette given 



by n^. We consider the case of a commensurate flux = 1/Q with Q integer. The 
magnetic held has opposite signs for opposite spin a ='[ (J,) = ±. This guarantees that the 
Hamiltonian is TRS invariant. Since = 1, the discrete translational invariance of the 
square lattice is broken down to translations by Q lattice sites on the y direction. We can use 
the previous notation for composite lattices with the relabeling iy) —)■ (ix, « + Qiy)- The 
Bloch functions and band dispersion are solutions of the Harper equation^^ (Supplementary 
Note 4). 

We are mainly interested in the limit of low hux density per plaquette = 1/Q 1. 

In this case, the bandwidth of each band is exponentially suppressed with respect to the 
band gap^^, thus Eq. (15) is satished. As shown in Supplementary Note 4, a self-consistent 
solution with Aq, = A = const. (A is now a real scalar) can be found, and therefore the 
result for the superfluid weight in Eq. (20) applies to the Harper-Hubbard model. The 
only missing piece is the evaluation of M. in Eq. (21)-(22). In the low magnetic field limit, 
a suitable approximation for the Bloch functions of the lowest bands consistent with the 
flat-band approximation ~ 

9n^.{a) ^ Y. , (25) 


where (pn{c() are the eigenfunctions of the harmonic oscillators if a is a continuous variable. 
In Supplementary Note 4, it is shown that, for the Harper model. 


M 


f2n + l —i \ 
y i 2n + lj ' 


(26) 


The superfluid weight in the h-th band (Landau level in the continuum) is proportional to 
oc 2n -|- 1. Note how the bound (23) is saturated for the lowest Landau level (all Landau 
levels in the continuum have^°’^^ \C\ = 1). By working directly in the continuum, we have 
obtained precisely the result contained in Eq. (20) and (26) (details are not provided here). 
More generally, Eqs. (20)-(21)-(22) are valid for a generic flat band, while Eq. (26) is specihc 
for the Harper Hamiltonian. 


III. DISCUSSION 

We have discovered that an invariant built from the quantum geometric tensor, which is 
intimately related to the Chern number, governs superfluidity in the flat-band limit. The 
inequality (23) implies that a topologically nontrivial flat band {C ^ 0) is guaranteed to 
have a hnite superfluid density in the presence of pairing in the system. Similar but more 
complicated bounds are also expected in three dimensions, since Aiij is positive semidehnite 
in general, and its imaginary part encodes three Chern numbers instead of one. This is the 
first time that the superfluid weight has been directly related to a topological invariant. 
Remarkably, BdG Hamiltonians with TRS and invariance under spin rotation around a 
given axis belong to the chiral unitary class AHI, whose ground state is topologically trivial 
in 2D according to the classification of Ref. 18, therefore we are referring to bulk superfluid 
transport and not to transport due to edge modes. 

In a flat band, mean-held theory is usually not adequate, however the BCS wavefunction, 
implicit in the Bogoliubov-de Gennes approach, is the exact ground state in the continuum 
limit of the Harper-Hubbard model considered here. This can be shown by mapping to the 
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wavefunction of a quantum Hall ferromagnet^*^^^® (see Methods IV D). Under this mapping, 
the result given by (20) and (26) for the superfluid weight of the Harper-Hubbard model 
translates into the spin stiffness or, equivalently, the counterflow-current superfluid density 
of a quantum Hall ferromagnet^^’^® with contact repulsive interactions. Whether mean-held 
theory can describe pairing in hat bands other than Landau levels is an open problem, 
analogous to the problem of characterizing the bands that can host a fractional Chern 
insulator^®’^^, but considerably less studied. We have checked that dynamical mean-held 
theory calculations (which treat local huctuations exactly) for the Harper-Hubbard model 
are indeed in excellent agreement with mean-held theory in the case of quasi-hat bands^^. 

Another problem of mean-held theory in 2D is that the transition to the normal state 
occurs at the Berezinsky-Kosterlitz-Thouless (BKT) transition temperature Tbkt, which is 
related to the superhuid density by a universal relation and is lower than the mean-held 
critical temperature. At half hlling, the estimated Tbkt is close to the mean-held transition 
temperature (see Supplementary Note 6 and Supplementary Figure 1) 

Tc = -Ar=o = > Tbkt • (27) 

Indeed, we hnd Tbkt ~ 0.25, 0.61, 0.75 Tc for h = 0,1, 2 respectively. 

The superhuid weight is a linear response transport coefficient, a ground state property, 
and it can be calculated exactly if the exact ground state is known^®, as in the case of the 
Harper-Hubbard model discussed above. As a consequence, it is not is necessary to employ 
beyond mean-held methods for estimating the superhuid weight^®. In summary, while the 
validity of mean-held theory for hat bands is in general an open question, the superhuid 
weight derived here for the Harper-Hubbard model is exact in the hat-band limit and a 
good approximation for quasi-hat bands. 

In ultracold gases, the atom-atom interaction is tunable, thus these systems are an ideal 
platform to conhrm our prediction that in a hat-band Dg oc U. In fact, it is possible to 
introduce complex hoppings in a lattice Hamiltonian (Peierls substitution) by Raman dress¬ 
ing^® or lattice shaking^®. Notably, the Harper model has been recently implemented with 
ultracold gases®®’^®. Whereas at the qualitative level superhuidity in ultracold gases is a 
well-established fact, a quantitative measurement of the superhuid weight has not been easy 
to perform so far. It has been proposed that the superhuid weight can be measured by an 
analog of the classic Adronikashvili experiments^, whereas the superhuid fraction of the uni¬ 
tary Fermi gas has been measured by means of second sound^^. Moreover, recent transport 
experiments with ultracold Fermi gases^^’S^ make it realistic to measure quantities like the 
superhuid weight. Currently, the main issue in ultracold gas experiments is the excessive 
heating present in experiments with artihcial gauge helds^®. Our estimates indicate (see 
Methods IV E) that superhuidity may be achieved in the future in topologically nontrivial 
hat bands that can be realized with ultracold atoms. Flat bands have been suggested as a 
possible mechanism to explain high-Tc superconductors'®’^^, and our results can be used to 
prove this hypothesis. If our results are generalized to the long-range Coulomb interaction, 
one more experimental context where they may be important are quantum Hall ferromagnets 
(c.f. the above-mentioned mapping of the superhuid weight [Eqs.(20) and (26)] to the spin 
stihness of a quantum Hall ferromagnet). In fact, a contact interaction is not an acceptable 
approximation in this case. 

Our results can be understood by distinguishing two possible ways to obtain a band 
of exactly degenerate states. On the one hand, the particles can be conhned in states 
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with negligible overlap by high potential barriers, or alternatively localization can occur 
in overlapping orbits due to (pseudo-)niagnetic helds or lattice geometry. In the latter 
case, the possibility of transport is a non-trivial question. The fact that we hnd a nonzero 
superfluid weight in a flat band can be understood by hnite overlap of the Cooper pairs, 
indeed pairing fluctuations support transport whenever Cooper pairs can be created and 
destroyed at distinct locations. Somewhat related in a work^^ that focused on condensation 
rather than superfluidity, an effective Hamiltonian for bosons in a flat band was derived 
by taking matrix elements of the interaction between overlapping Wannier functions, which 
produced an effective hopping for the particles. In the work of Provost and Vallee^®, pointing 
out for the hrst time the natural geometric structure present in a manifold of quantum states, 
it is suggested that macroscopic quantum systems that exhibit collective behaviour might 
be those where the quantum metric has direct physical signihcance, an intuition that has, 
in some sense, materialized in our results which showed the connection between quantum 
metric and superfluid weight. It is an intriguing topic for future research to understand 
whether the pairing fluctuations and macroscopic phase of a superfluid have any connection 
to the fact that the quantum metric equals the fluctuations in the quantity that generates 
the path of a quantum state in the manifold^®. 

While the above discussion may help to guide the intuition, the rigorous framework for 
future work is given by our results on the important role of Wannier functions in superfluid 
transport. As we have shown, the bound on the superfluid weight translates into a bound 
on the the localization functional for Wannier functions^^. A nonzero Chern number implies 
that the Wannier functions have algebraically decaying tails^^, and this explains the bound 
Dg > \C\. But the Wannier functions can also be delocalized on a short range only, which 
is consistent with the fact that the superfluid weight is related to an invariant distinct from 
the Chern number. In general, we propose (quasi-)flat bands as a viable way to increase the 
critical temperature in novel superconducting materials, while at the same time preserving 
the dehning properties of superconductors. We expect the invariant M. that controls the 
superfluid weight in a flat band to play a central role in this research effort. 


IV. METHODS 

A. Derivation of the Bogoliubov-de Gennes Hamiltonian 

The BdG Hamiltonian in Eq. (6) is important for our purposes, and here we clarify 
its derivation. The hopping matrix has the same discrete translational symmetry as the 
Bravais lattice, since j). By expanding the held operators into plane waves 

c\aa = Xlk the kiuetic term of the Hamiltonian can be block-diagonalized in 

momentum space 


(k)]Q.,/3C^ko- 

(28) 

i«J/3 

k,a,/5 


with 


(29) 


i-j 


It is convenient to introduce a Nambu spinor Ck = built out of the operators 

Caku in the plane wave basis (see Eq. (28)) and write the mean-held Hamiltonian as "Hm.f. = 
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Ek4^k(q)ck with 


ir^(k-q) -/il 

A 


(30) 


^k(q) = 


-(A't(k + q)-Al)) ^ 


To cast the mean-held Hamiltonian in Nambu form, we have anticommuted the spin-down 
creation and annihilation operators in the kinetic energy term and used TRS in the form 
(A'*'(k)) = K^{—'k). All the c-number terms in the mean-held Hamiltonian have been 

dropped since they do not ahect the superhuid weight (see Supplementary Notes 1 and 3). 
A further canonical transformation is performed to go from the basis given by the orbitals 
within a unit cell, labelled by a, j3, to the basis that diagonalizes the kinetic Hamiltonian, 
that is, the Bloch functions labelled by n. More precisely, the transformation reads c„kt = 

E„[^{k-q)t]a,nt^nkt and = En[^(-k-qU]^,„dEk4 = En[%+q)t]/3,n4- 14 - t^is Way, 

Eq. (6) is obtained. 


B. Definition of a generic function of an Hermitian matrix 

In Eq. (10) and Eq. (11), the absolute value | ■ | and the sign function sign(-) of the BdG 
Hamiltonian hfk(q) are used. In general, a function /(■) of an Hermitian matrix H = UDU\ 
diagonalized by the unitary matrix U and by the real diagonal matrix D, is dehned as the 
function of the eigenvalues f{H) = Uf(D)W. 


C. Positive semidefiniteness of the quantum geometric tensor 


In Eq. (22), the projection Q\^ of the unitary matrix Q\^ on the n-th band is dehned by 


[^k]o,l — fi'nk(ci) — [^k] 


Q;,n ? 


^k^k = 1 , 


= PfL 


nk • 


(31) 


Pfik is a projection operator, a positive semidehnite {Pnk > 0) and idempotent {P?i^ = Pnk) 
operator. The matrix Qk is just a column vector in Eq. (31), but it can be a rectangular 
matrix for a group of degenerate hat bands, for example. Since the dispersion is hat, 
^k characterizes the hat band completely. The positivite semidehniteness of the projector 
Pnk and of its complement 1 — P„k implies that the matrix Hjj(k) in Eq. (22) is positive 
semidehnite since it can be written in the form 


%(k) 


2Tr 


(4.^i)(i 


Pnk){,dkjQk) 


(32) 


The invariant matrix in Eq. (21) is also positive semidehnite Ad > 0 since it is a linear 
combination with positive coefficients of the positive semidehnite matrices i3jj(k). Interest¬ 
ingly, the Berry curvature (Imi3jj(k)) and the Chern number of a set of bands are obtained 
by adding the respective contributions of all bands in the set, whereas the quantum metric 
(Rei3jj(k)) is not additive due to the second term in Eq. (22), which is real and involves a 
double sum over the band index. 


D. Exactness of the BCS wavefunction 

The BCS wavefunction can be shown to be the exact ground state of the Harper-Hubbard 
model in the hat-band limit. To take the hat-band limit of the Harper-Hubbard model, it is 
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necessary to take the limit of low magnetic flux. The problem is mapped into that of particles 
in the continuum in the presence of a constant magnetic flux (Landau problem). We consider 
a general form for the interparticle interaction potential ld(r) = (27r)~^ f d^qv(q)e^^'’^ and 
perform the projection of the interaction term into the h-th Landau leveh^ 

^ ^ _ ( 33 ) 

Here £b is the magnetic length, Ln{x) is the n-th Laguerre polynomial and p„qo- = p|j(_q)o- 
are projected density operators that obey the Girvin-MacDonald-Platzman algebra 



where p A q = pxQy — PyQx- In the Landau gauge, the explicit expression for the projected 
density operators is 

, (35) 

k 

The annihilation (creation) operators are labelled by the Landau level index n and the 
momentum k along the x direction which is conserved in the Landau gauge. Notice that if 
n(q) > 0, the interaction Hamiltonian (33) is repulsive between particles with the parallel 
spins and attractive between particles with antiparallel spins. It is straightforward to verify 
that the operator p^^q^ — in the Hamiltonian annihilates the BCS wavefunction 

|BCS) = JJ(m + |0 ) (36) 

k 

for arbitrary values of u and v and any value of q, that is, the BCS wavefunction is a zero 
eigenvector of the Hamiltonian. Normalization requires that \u\^ + |np = 1. A possible 
parametrization is m and v = e^^\/l — v with v the hlling and and arbitrary 

phase. Since the Hamiltonian (33) is a positive semidehnite operator for t;(q) > 0, the BCS 
wavefunction must be the ground state since it is a zero eigenvector. 

An alternative way to interpret this result is well known in the context of quantum Hall 
physics^^’^®. By performing a particle-hole transformation of the form c^, the BCS 

wavefunction is transformed into the wavefunction of a completely polarized ferromagnet 

iFerro) = JJ (ucl^ + vcl^ |0) . (37) 

k 

This is a simple Slater determinant where all the states with spin wavefunction n |j') -f- 
u |j,) are occupied. Under the same transformation, the interparticle interaction becomes a 
repulsive interaction, which is completely isotropic in spin space. It is easy to understand 
why the wavefunction Eq. (37) is the ground state. According to Hund’s rule, the interaction 
energy is minimized if the all the spins are parallel (a consequence of the Pauli exclusion 
principle), and in a Landau level, there is no kinetic energy cost that prevents a complete 
alignment. Indeed, this extreme ferromagnetic state has been observed in experiments in the 
quantum Hall regime^"^. It is important to note that the 2 : component of the magnetization 
in the ferromagnetic state is mapped by the particle-hole transformation into the total 
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number of particles on the superconducting side (and vice versa). Therefore, whereas the 
wavefunction (37) is the ground state when a spinful Landau level is half-hlled, the BCS 
wavefunction is the correct ground state for any hlling. 

In the limit of a contact interaction, the repulsive interaction between particles with 
parallel spins disappears and one is left with a purely attractive interaction, that is, the 
continuum limit of the Harper-Hubbard model considered here. 


E. Estimate of the critical temperature for the Harper-Hubbard model 

To estimate the critical temperature for an actual ultracold gas experiment, we consider 
fermionic ®Li atoms in an optical lattice with a typical wavelength of the laser standing wave 
A = 1064 nm = 2T[/k and the corresponding recoil energy given by E^. = {hkY/2m6i^i = 
1.4 /iK. The hopping energy scale J can then be estimated from the approximate formula 



Using the same ratio Vq/E^ ^ 7 as in Ref. 29 between the amplitude Vq of the optical 
lattice potential and the recoil energy, one obtains J ~ 70 uK. In Supplementary Note 
4 and Supplementary Figure 2, we estimate that in the isolated flat-band approximation 
for the time-reversal invariant attractive Harper-Hubbard model, the mean-held critical 
temperature is of the order of k-^T^ ~ 0.02 J, which implies a BKT transition temperature in 
the order of the nanoKelvin. Such a low temperature results just because we wished to be 
able to use the analytical results derived here, which requires pairing within a single band and 
thus U needs to be smaller than the gaps to neighboring bands, Eq.(15). Conceptually the 
same results can, however, be achieved when several (but not all) hat (or nearly hat) bands 
of the composite bands participate in pairing, only that the theoretical analysis becomes 
more involved. Then, the limit on U is relaxed, and Tc can be substantially increased. 


* paivi.torma@aalto.fi 

^ Kopnin, N. B., Heikkila, T. T. &: Volovik, G. E., High-temperature surface superconductivity 
in topological flat-band systems, Phys. Rev. B 83, 220503(R) (2011). 

^ Khodel, V. A., &: Shaginyan, V. R., Superfluidity in systems with fermion condensate, JETP 
Letters 51, Issue 9, 553-555 (1990) [Pis’ma Zh. Eksp. Teor. Fiz. 51, No. 9, 488-490 (1990)]. 

^ Khodel, V. A., Shaginyan, V. R. Khodel, V. V., New approach in the microscopic Fermi 
systems theory, Phys. Rep. 294, 1-134 (1994). 

^ Volovik, G. E., A new class of normal Fermi liquids, JETP Letters 53, Issue 4, 222-225 (1991) 
[Pis’ma Zh. Eksp. Teor. Fiz. 53, No. 4, 208-211 (1991)]. 

^ Heikkila, T. T., Kopnin, N. B. Sz Volovik, G. E., Flat bands in topological media, JETP Letters 
94, Issue 3, 233-239 (2011). 

® Nandkishore, R., Levitov, L. S. &: Chubukov, A. V., Ghiral superconductivity from repulsive 
interactions in doped graphene. Nature Phys. 8, 158163 (2012). 

^ Kopnin, N. B., Ijas, M., Harju, A. & Heikkila, T. T., High-temperature surface superconduc¬ 
tivity in rhombohedral graphite, Phys. Rev. B 87, 140503(R) (2013). 


14 



® Uchoa B. &: Barlas, Y., Superconducting states in pseudo-Landau-levels of strained graphene, 
Phys. Rev. Lett. Ill, 046604 (2013). 

® Tang, E., &: Fu, L., Strain-induced partially flat band, helical snake states and interface super¬ 
conductivity in topological crystalline insulators, Nature Phys. 10, 964969 (2014). 

Yudin, D., Hirschmeier, D., Hafermann, H., Eriksson, O., Lichtenstein, A. 1. &: Katsnelson, M. 
L, Fermi condensation near van Hove singularities within the Hubbard model on the triangular 
lattice, Phys. Rev. Lett. 112, 070403 (2014). 

Khodel, V. A., Clark, J. W., Popov, K. G. &: Shaginyan, V. R., Occurrence of flat bands in 
strongly correlated Fermi systems and high-Tc superconductivity of electron-doped compounds, 
JETP Letters 101, Issue 6, 413-418 (2015). 

Scalapino, D. J., White, S. R. & Zhang, S.-C., Superfluid density and the Drude weight of the 
Hubbard model, Phys. Rev. Lett. 68, 2830-2833 (1992). 

Scalapino, D. J., White, S. R. &: Zhang, S.-C., Insulator, metal, or superconductor; The criteria, 
Phys. Rev. B 47, 7995-8007 (1993). 

Grosso, G. &: Parravicini, G. P., Solid State Physics, 2nd ed., Elsevier (2014). 

Wannier, G. H., Dynamics of band electrons in electric and magnetic fields. Rev. Mod. Phys. 
34, 645-655 (1962). 

Nenciu, G., Dynamics of band electrons in electric and magnetic fields; rigorous justification of 
the effective Hamiltonians, Rev. Mod. Phys. 63, 91-127 (1991). 

Kopnin, N. B., Surface superconductivity in multilayered rhombohedral graphene; supercurrent, 
JETP Letters 94, Issue 1, 81-85, (2011). 

Schnyder, A. P., Ryu, S., Furusaki, A. & Ludwig, A. W. W., Classification of topological 
insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008). 
Kane, C. L. &: Mele, E. J., Z 2 Topological order and the quantum spin hall effect, Phys. Rev. 
Lett. 95, 146802 (2005). 

Thouless, D. J., Kohmoto, M., Nightingale, M. P. & den Nijs, M., Quantized Hall conductance 
in a two-dimensional periodic potential, Phys. Rev. Lett. 49, 405-408 (1982). 

Hasan, M. Z. & Kane, C. L., Colloquium; Topological insulators. Rev. Mod. Phys. 82, 3045-3067 

( 2010 ). 

Qi, X.-L. (fe Zhang, S.-C., Topological insulators and superconductors. Rev. Mod. Phys. 83, 
1057-1110 (2011). 

Bernevig, B. A. with Hughes, T. L., Topological Insulators and Topological Superconductors, 
Princeton University Press (2013). 

Marzari, N., Mostofi, A. A., Yates, J. R., Souza, 1. & Vanderbilt, D., Maximally localized 
Wannier functions; Theory and applications. Rev. Mod. Phys. 84, 1419-1475 (2012). 

Brouder, C., Panati, G., Calandra, M., Mourougane, C. & Marzari, N., Exponential Localization 
of Wannier Functions in Insulators, Phys. Rev. Lett. 98, 046402 (2007). 

26 Provost, J. P. & Vallee, G., Riemannian structure on manifolds of quantum states, Commun. 
Math. Phys. 76, 289 (1980). 

Berry, M. V., “The quantum phase, five years after”, pp. 7-28 in Geometric phases in physics, 
eds. A. Shapere & F. Wilczeck , World Scientific (1989). 

Wang, L., Hung, H.-H. &: Troyer, M., Topological phase transition in the Hofstadter-Hubbard 
model, Phys. Rev. B 90, 205111 (2014). 

Aidelsburger, M., Atala, M., Lohse, M., Barreiro, J. T., Paredes, B. & Bloch, L, Realization 
of the Hofstadter Hamiltonian with ultracold atoms in optical lattices, Phys. Rev. Lett. Ill, 
185301 (2013). 


15 



Miyake, H., Siviloglou, G. A., Kennedy, C. J., Burton, W. C. & W. Ketterle, Realizing the 
Harper Hamiltonian with laser-assisted tunneling in optical lattices, Phys. Rev. Lett. Ill, 
185302 (2013). 

Taylor, E., Griffin, A., Fukushima, N. & Ohashi, Y., Pairing fluctuations and the superfluid 
density through the BGS-BEC crossover, Phys. Rev. A 74, 063626 (2006). 

Suhl, H., Matthias, B. T. &: Walker, L. R., Bardeen-Cooper-SchriefFer theory of superconduc¬ 
tivity in the case of overlapping bands, Phys. Rev. Lett. 3, 552-554 (1959). 

Nagamatsu, J., Nakagawa, N., Muranaka, T., Zenitani, Y. & Akimitsu, J., Superconductivity 
at 39K in magnesium diboride. Nature 410, 63 (2001). 

Xi, X. X., Two-band superconductor magnesium diboride. Rep. Prog. Phys. 71, 116501 (2008). 
Chang M.-C. & Niu, Q., Berry phase, hyperorbits, and the Hofstadter spectrum, Phys. Rev. 
Lett. 75, 1348-1351 (1995). 

Chang, M.-C. & Niu, Q., Berry phase, hyperorbits, and the Hofstadter spectrum: Semiclassical 
dynamics in magnetic Bloch bands, Phys. Rev. B 53, 7010-7023 (1996). 

Kopnin, N. B. & Sonin, E. B., BCS Superconductivity of Dirac Electrons in Graphene Layers, 
Phys. Rev. Lett. 100, 246808 (2008). 

Kopnin, N. B. &: Sonin, E. B., Supercurrent in superconducting graphene, Phys. Rev. B 82, 
014516 (2010). 

Neupert, T., Chamon C. &: Mudry, C., Measuring the quantum geometry of Bloch bands with 
current noise, Phys. Rev. B 87, 245103 (2013). 

Dobardzic, E., Milovanovic, M. V. &: Regnault, N., Geometrical description of fractional Chern 
insulators based on static structure factor calculations, Phys. Rev. B 88, 115117 (2013). 

Roy, R., Band geometry of fractional topological insulators, Phys. Rev. B 90, 165139 (2014). 
Marzari, N. &: Vanderbilt, D., Maximally localized generalized Wannier functions for composite 
energy bands, Phys. Rev. B 56, 12847 (1997). 

Harper, F., Simon, S. H. &: Roy, R., Perturbative approach to flat Chern bands in the Hofstadter 
model, Phys. Rev. B 90, 075104 (2014). 

Girvin, S. M. &: MacDonald, A. H., “Multicomponent Quantum Hall Systems: The Sum of 
Their Parts and More”, p. 161 in Perspectives in Quantum Hall Effects, eds. S. das Sarma Sz 

A. Pinczuk, Wiley-VCH (1996). 

Spielman, 1. B., Eisenstein, J. P., Pfeiffer, L. N. & West, K. W., Resonantly enhanced tunneling 
in a double layer quantum Hall ferromagnet, Phys. Rev. Lett. 84, 5808-5811 (2000). 

Joglekar, Y. N. &: MacDonald, A. H., Microscopic functional integral theory of quantum fluc¬ 
tuations in double-layer quantum Hall ferromagnets, Phys. Rev. B 64, 155315 (2001). 

Vanhala, T.L, Peotta, S., Troyer, M. & Torma, P., in preparation. 

Evans, W. A. B. & Rashid, R. 1. M. A., A conserving approximation evaluation of superfluid 
density within the pair theory of superfluids, J. Low Temp. Phys. 11, 93-115 (1973). 
Jimenez-Garcfa, K., LeBlanc, L. J., Williams, R. A., Beeler, M. C., Perry, A. R. & Spielman, 1. 

B. , Peierls substitution in an engineered lattice potential, Phys. Rev. Lett. 108, 225303 (2012). 
Struck, J., Qlschlager, C., Weinberg, M. , Hauke, P., Simonet, J., Eckardt, A., Lewenstein, M., 
Sengstock, K. & Windpassinger, P., Tunable gauge potential for neutral and spinless particles 
in driven optical lattices, Phys. Rev. Lett. 108, 225304 (2012). 

Cooper, N. R. &: Hadzibabic, Z., Measuring the superfluid fraction of an ultracold atomic gas, 
Phys. Rev. Lett. 104, 030401 (2010). 

Sidorenkov, L. A., Tey, M. K., Grimm, R., Hou, Y.-H., Pitaevskii, L. &: Stringari, S., Second 
sound and the superfluid fraction in a Fermi gas with resonant interactions. Nature 498, 78-81 


16 



(2013). 

Stadler, D., Krinner, S. , Meineke, J., Brantut, J. P. &: Esslinger, T., Observing the drop of 
resistance in the flow of a superfluid Fermi gas, Nature 491, 736739 (2012). 

Krinner, S., Stadler, D., Husmann, D., Brantut, J. P. & Esslinger, T., Observation of quantized 
conductance in neutral matter. Nature 517, 6467 (2015). 

Huber, S. D. & Altman, E., Bose condensation in flat bands, Phys. Rev. B 82, 184502 (2010). 

Acknowledgments 

This work was supported by the Academy of Finland through its Centres of Excellence Pro¬ 
gramme (2012-2017) and under Project Nos. 263347, 251748, and 272490, and by the Euro¬ 
pean Research Council (ERC-2013-AdG-340748-CODE). We acknowledge useful discussions 
with Jami Kinnunen, Jildou Baarsma and Tuomas Vanhala. We thank Antti Paraoanu for 
producing the illustrations for hgures la)-lc) and 2c)-2d). 

Authors Contributions 

All authors contributed to all aspects of this work. 

Competing financial interests 

The authors declare no competing hnancial interests. 


17 




d) 

A^ij — + icijC ^ 0 

\ I 

D, > ICI 

FIG. 1. (a) Localized Wannier functions are obtained from the Bloch functions of a set of bands, 

called a composite band^^’^^. To have superfluidity in a flat band, the pairing takes place only in 
a subset of the bands within the composite band, e.g. in a single flat band. While the Wannier 
functions built from the Bloch functions of the band where pairing takes place are delocalized 
due to the nonzero Chern number^^ C ^ 0 (b) the Wannier functions of the composite band are 
exponentially localized (c). We show that the superfluid weight Dg in a flat band is given by the 
Brillouin-zone average of the quantum metric^®’^^, which is the real part of an invariant Aiij. The 
imaginary part of Aiij gives the Chern number C. The positivite semidefiniteness of Aiij leads to 
the bound Dg > IC*! on the superfluid weight. 
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FIG. 2. (a) We consider lattices in three and two dimensions (2D in the figure). Periodic bound¬ 
ary conditions (PBC) are used. The lattice contains W unit cells and each unit cell contains A^orb 
sites/orbitals. The vectors aj {i = 1, 2,3) are the fundamental vectors of the Bravais lattice^^ while 
the vectors (a = 1,..., Norh) ure the positions of the centers of the orbitals (Wannier functions) 
within a unit cell. A single orbital is specified by a triplet (or pair) of integers i = {ix, iy, iz) and by 
the sublattice index a and is centered at the position vector rio = ixai + iya 2 + iza^ + be,, (b) The 
band structure is obtained by solving the Schrodinger equation + f^(r)^ '4’nk = 

with periodic potential P(r) = I/(r -|- aj). It consists of the band dispersions e„k) with n the band 
index and Kk the lattice quasimomentum, and the periodic Bloch functions (7nk(r) = 5'nk(r + Uj) 
(Bloch functions for brevity) obtained from the Bloch plane waves '0nk(r) = e**^'“'5„k(r)- We con¬ 
sider a composite band, that is, a subset S of contiguous bands well separated in energy from other 
bands. The Chern numbers Cn for individual bands calculated from the Bloch functions may be 
nonzero (such as the flat band n = 2 in the figure), but their sum equals zero X]ne5 ~ 

Chern number refers to spin-resolved bands since the spin along a quantization axis (conventionally 
the z axis) is conserved, (c) The Wannier functions, defined as the Fourier transform of the Bloch 
functions, allow us to derive a tight-binding Hamiltonian that reproduces exactly a single band or a 
composite band of the original continuum Hamiltonian (see Supplementary Note 2). Since individ¬ 
ual bands may be topologically nontrivial with nonzero Chern numbers, their Wannier functions 
Wni'r) are not exponentially localized^®, and the Peierls substitution in the effective Hamiltonian 
is therefore not justified^®’^®. (d) By constructing Wannier functions as linear superpositions of 
Bloch waves of all bands in the composite band, exponentially localized Wannier functions rcQ(r) 
can be created. The mixing of the different bands is provided by the unitary matrix Ua,n(M)- This 
justifies the Peierls substitution for a composite band S. 
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Supplementary Figures 
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Supplementary Figure 1: On the left plot the energy gap A(r)/A(0) (blue solid line) and the stiffness 
Tra(T)/A(0) (red dashed line) as a function of temperature are show. The stiffness Tn is shown for several 
values of h. The intercepts between the black straight line and the red dashed lines correspond to the 
solutions of Eq. (95) which gives the BKT temperature Tbkt for several values of n. The calculated 
solutions of Eq. (95) are shown in the right plot. Only half-filled flat bands {v = 1/2) are considered here. 
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Supplementary Figure 2: In the upper panels the band structure of the Harper model for three different 
values of the number of magnetic flux quanta per plaquette = 1/Q = gj gj | is shown. The band 
dispersion is shown as a function of the wavevector kx for twenty different values of ky. Both kx, ky are 
restricted to the reduced Brillouin zone kx, ky e [—tt/(Q a), tt/ {Qa)] due to the degeneracy of the dispersion 
in the Harper model (53). Notice how for decreasing the lower bands approach the flat-band limit. In 
the lower panel the corresponding energy gap A calculated self-consistently by using Eq. (59)-(60) (solid 
lines) is shown as a function of the normalized total particle density rip/ (2n^) = n + u. The value of the 
coupling constant U/ J in all three cases is approximately equal to one-fourth the energy gap between the 
two lowest bands, more precisely U/J = 0.31 for = 1/3, U/J = 0.39 for = 1/5 and U/J = 0.35 for 
— 1/1. The critical temperature k^T^ = = At=o/ 2 shown in the lower panels is the mean-field 

critical temperature in the isolated flat-band limit. The numerical self-consistent solution of Eq. (59)-(60) 
is compared to the approximate solution in the isolated flat-band limit (Eq. (61)-(64)) (dashed lines) for 
three different temperatures T/T^ = 0.2,0.75,0.9. The isolated flat-band approximation is not good for 
any band in the case = 1/3, but it can be used for the first band for n^f, = 1/5 and for the first two 
bands for = 1/7. 
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Supplementary Note 1: Superfluid weight from the grand potential 

The aim of this Supplementary Note is to justify Eq. (1) in the main text stating that the superfluid weight 
can be calculated by taking successive derivatives of the grand potential A, q) only with respect to 
the wavevector q. The gap function A(q) and the chemical potential /i(q) are themselves a function of q 
(for fixed particle density rip). This dependence is in principle important and requires the solution of the 
self-consistency equations in A and /r for nonzero values of q. It is shown here that this not the case. In 
the following we do not make use of translational invariance since the result that we want to prove holds 
in general. The greek indices a,j3 that appear in the case of a composite lattice will be suppressed as well 
since what we are going to prove holds equally well in single-band and multiband systems. The mean-held 
decoupling of the Hubbard interaction is (including the c-number term) 

j j 

It will be shown below that the self-consistency equation for the pairing order parameter (also called energy 
gap function) Aj = —U{cjiCj-f) can be obtained as an extremizing condition for the grand potential. 

The mean-held Hamiltonian T-Lm.f., comprising the single particle (quadratic) term in Eq. (3) in the 
main text and the mean held approximation (1) for the interaction term can be conveniently rewritten in 
matrix form in Nambu space 

= c'l'HBdG(q)c-f TrH'^(q) -b ^ ^ |Ai(q)|2 , (2) 

where the following dehnitions have been used 



(Nambu spinor), 


-ffBdG(q) 


/ Hj-iq) A;(q)(Iij\ 

VA?(q)(5;j -Hi^j(q)y ' 


^fij(q) = A:^(q) - M(q)<^ij + , 

= A:Jj(-q) - M(q)<JiJ + , 


( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 


The trace in Eq. (2) comes from the constant term produced when anticommuting fermion operators 
~ ('^iJ Kronecker delta function). In our case we have anticommuted the fermionic 

operators relative to the down spin. 

In the above equations several quantities depend on the wavevector q. The hopping matrix itrfj(q) = 
depends on q according to the Peierls substitution, Eq. (2) in the main text. The pairing 
potential A;(q) is a function of q due to the fact that the mean-held Hamiltonian (2) must be solved 
self-consistently, while the chemical potential depends on q since the particle number per lattice site rip 
is hxed. We have introduced an additional spin-independent scalar potential V) in the diagonal terms of 
^^BdG (see Eq. (5)-(6)) since our derivation applies to a system without translational invariance. In Eq. (6) 
we have used the time-reversal symmetry (TRS) of the hopping term (it'j4j(q))* = it'4^(—q). Notice that a 
hnite wavevector q breaks TRS since it has the effect of inducing a hnite supercurrent in the ground state. 
However TRS connects the Hamiltonians with opposite wavectors q —)■ —q, in fact under TRS the sign of 
the supercurrent is reversed. 

A basic thermodynamic function we are interested in is the thermodynamic grand potential dehned 
from the grand partition function Zq 


n{T, /.(q), Ai(q), q) = -1 InZj, = -^ lnTY[e-^«-* ] = TYH;(q) + ^ 


|Ai(q)P 

U 


— — lnTr[e 


-/3ctifBdG(q)c] 


(7) 


3 



The grand potential is a function of the wavevector q either directly, since q appears in the diagonal 
blocks of i?BdG(Q) through the kinetic energy term itr'^(q) [see Eq. (4)] or indirectly through the quantities 
/^(q), Ai(q). This distinction is important for what follows. The last term on the right hand side of the 
above equation can be evaluated by diagonalizing the BdG Hamiltonian. We use the notation 

/fBdG(q) = W(q)E(q)Wt(q). (8) 

The matrix W(q) is unitary, while E{q) is diagonal and the diagonal elements are the eigenvalues ^^(q), 
namely 

[-E(q)]a,b = Ea{q}da,b = £'a(/^(q), Ai(q),q)(5a,fc. (9) 

The indices a,b = 1,, 2Ns run over the 2Ns eigenvalues of the BdG Hamiltonian and Ng = ApAorb is 
the number of lattice sites (see Fig. 2 in the main text). The square bracket notation [M]a,b denotes the 
matrix element in row a and column b of the matrix M while 5afi is the Kronecker delta. The eigenvalues 
£'a(q) of the BdG Hamiltonian depend on q either directly or indirectly in the same way as the grand 
potential. The notation Ea{<\) is just a convenient short-hand for the right hand side of Eq. (9). 

The canonical transformation of the fermionic operators c 


c = Vy(q)7, (10) 

preserves the anticommutation relations. The operators 7 = ( 71 ,... ,j 2 Ns)’^ are therefore annihilation or 
creation operators of the fermionic quasiparticles that are the elementary excitations above the ground 
state of the mean-field Hamiltonian. Using the canonical transformation (10) the trace on the right hand 
side of Eq. (7) can be evaluated 

rpj.jg-/3ctHBdG(q)e] ^ = TT Tr[e“^^“("‘)^“^“] = 17(1 -F . (11) 


The above identity is a consequence of the fact that the number operators ylya are mutually commuting 
and that the trace of a tensor product is the product of traces Tr[H 0 i?] = Tr[yl]Tr[i?]. The grand potential 
thus reads 

0(r,M(q), A;(q),q) = Tr A;(q) + ^ . ( 12 ) 

i ^ a 

At equilibrium the grand potential attains a minimum, therefore it must be stationary with respect to 
variations of the pairing potential, namely 

^ ^ ^ ^ ^ Tr[pthC 4 C:t] = ^ + (citCit). (13) 


The expectation values {•) = Tr[pth • ] are taken with respect to the equilibrium state described by the 
thermal density matrix pth = f Eq. (13) is the equation that guarantees the self-consistency of 

the mean-field approximation and is written in the most general form valid for an inhomogeneous system. 
The total particle number A can be fixed using the thermodynamic relation 


A = - 


dn 

dfi 




(14) 


Eqs. (13)-(14) must in general be solved simultaneously by numerical means. Suppose that such a solution 
has been found for a given q. We now want to show that a self-consistent solution of the time-reversal 
transformed state with q —> — q can be immediately obtained according to the transformation 


Ai(q) = A?(-q), 

M(q) = M(-q) • 
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(15) 

(16) 



The proof proceeds by first solving the BdG Hamiltonian -HBdG(“Ct) assuming the above relations and 
then showing that indeed the self-consistency equations (13) and (14) are satisfied. Eqs. (15)-(16) imply 
that the BdG Hamiltonian has the following symmetry property 

tyHBdG{^)ty = -HBdG(-q) with ty = 0 Iat, • (17) 

The symmetry operator ty is simply a Pauli matrix that acts on the particle-hole space of the BdG 
Hamiltonian, namely the matrix blocks in Eq. (4), and acts as the identity on the orbital degrees of 
freedom (lattice sites). Eor q = 0 this symmetry has the properties 

= tytl = l (18) 

and is called a chiral or sublattice symmetry according to the terminology of Ref. 1, but it is not related to 
a particular transformation of the lattice in our case. The symmetry (17) is the result of the particle-hole 
symmetry which is always present in a BdG Hamiltonian together with time-reversal symmetry (TRS) and 
invariance with respect to rotations of the spin in the plane perpendicular to the z axis. TRS is broken for 
q 7 ^ 0 as discussed before. This symmetry implies that for q = 0 our system belongs to the AHl (chiral 
unitary) symmetry class, while for q 0 the system belongs to the A (unitary) symmetry class^. 

The symmetry in Eq. (17) implies that 


tyW{q_)ty = W(-q) and tyE{-q)ty = -E{q) . 


(19) 


It is important to note that the diagonalization of the BdG Hamiltonian i7BdG(q) is not uniquely defined 
and we can choose to permute the eigenvalues according to the transformation E{q) —> PE{q)P and 
W(q) —> W(q)P with P a permutation matrix. However if we require the matrix elements of W(q) and 
E(q) to be continuous and differentiable functions of q then the choice (19) is the only possible one. 

In the following we make use of the following additional matrix notations. A diagonal matrix D — 
diag((ia) is completely characterized by the set of its diagonal matrix elements {da}a=i,...,dim(D)- A generic 
function / of a diagonal matrix D = dmg{da) is taken element by element f{D) = diag(/(da)) and a 
function of a Hermitian matrix H = UDU^, with U the unitary matrix that diagonalizes H, is defined by 
f{H) = Uf{D)U^. In particular we use the modulus of an Hermitian matrix |i7| = U\D\U''^ and the sign 
function sign(i7) = U sign(Il) t/i. The sign function is defined in a such a way that sign(i7)|i7| = El. 

We can use the above identities (19) to calculate the expectation values of bilinear combinations of the 
fermionic operators Cio.,Cj^ which can be conveniently organized into a matrix 


C 0 C )_ 


t\ - 


= yy(-q)(7 0 7'f)_qyy^(-q) = w(-q)^p^g^p^y^vy^(-q) 


= 1 


AitCj.|.)-q (CitCj4,)-q 

- f^W(q) w^'(q)t^ = 1 - ty{c 0 ct)qt 


(Gt%)q 


( 20 ) 


Since the quasiparticle operators % diagonalize the BdG Hamiltonian Et^dGi^) the matrix (7 0 7 ^)q is 
diagonal 

(7® 7*), = diag + + 1 ■ 

Specializing the result in Eq. (20) for i = j one obtains Eqs. (15). Moreover since the total particle 
number N — unchanged also the chemical potential ^(q) is an even 

function, validating Eq. (16). This means that from a self-consistent solution for a certain value of q we 
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can immediately provide a self-consistent solution for the TRS conjugate ground state q —> —q. Eq. (20) 
shows that in principle it is possible to have a nonzero spin density in a state with finite current. From 
Eq. (20) the density of the z component of the spin S'^(q) = ni'i'(q) — 'Ui|(q) = (4fCi'|')q — (c;|.Cit)q is an 
odd function of q, which is consistent with TRS. 

The properties (15)-(16), imply that the grand potential is an even function of q. In fact, consider only 
the first and the last term in Eq. (12), since the remaining term is an even function of q according to our 
previous discussion, 

Tt ^ ln(l -|- = Tr H^{q) — ^ ln(l -|- 

a (22) 

= Trgi'(q) - TrgBdG(q) " ^ ^^1^(1 -|- e = Tr iji'(q) - - ^ ln(l-|-e . 


This implies that q = 0 is a stationary point of the grand potential, a candidate ground state. In order to 
check the stability of this stationary point it is necessary to calculate the second derivative with respect 
to q. This amounts to the evaluation of the superfluid weight Ds- The superfluid weight is defined as 
the second derivative of the free energy E(r, iV, A;(q), q) = n(T,/x(q), Ai(q), q) -|- A^^(q) since the total 
particle number is constant. 


. , ^ 1 d‘^F ^ 1 d‘^{n + piN) 

® dqidqj Vh? dqidqj 


(23) 


where V is the system volume (in two dimensions the area A is used instead). In general the superfluid 
weight in an anisotropic system depends on the direction and is thus a tensor. As emphasised in our 
notation the quasiparticle energies and thus the grand potential depend on q both through /i(q) and A(q) 
and directly through the diagonal terms of the BdG Hamiltonian. It will be now shown that in fact it 
is not necessary to know the dependence of /x(q), A(q) to calculate the superfluid density and only their 
q = 0 value is needed. This result has been proved in Ref. 2, but we repeat the proof here with some 
modifications. By a suitable gauge transformation Cjo- —)• Cia- it is possible to take Ai(q) real and at 

the same time preserving the property (15). This gauge will be used in the following. Taking the total 
derivative of the free energy with respect to q gives 


OF 

dqi 


E 


Mi 



< dfi 

-h 1 

dfj, 

77^ + 

dn 


dn 



A.q ) 

dqi 

dqi 

A,./x 

dqi 


(24) 


This result is valid for arbitrary q and it relies on the fact that p. and A; satisfy Eqs. (13)-(14), i.e. they 
are self-consistent solutions. As a consequence a stationary point of the free energy corresponds to the 
condition dil/dqi\^, ^ = 0 on the grand potential which is satisfied for q = 0 as we have seen above. From 
Eqs. (15)-(16) and the fact that Ai(q) is real, one has dAi/dqi = dp,/dqi = 0 at q = 0. This last result is 
useful when taking one more derivative of the free energy and evaluating it at equilibrium 


d‘^F 

d^Q dAi 

1 d'^fl dfi 

dfl 

1 

dn 


dqjdqi 

q=o i 9Aidqi dqj 

q =0 9ndqi dqj 

q=0 9qjdq^ 

q^o dqj dqi 

q=0 


In the above equation the partial derivatives of the grand potential with respect to A;, /i, qi are always taken 
keeping the other quantities fixed (introducing a notation to specify this would be quite cumbersome). The 
derivatives of the free energy with respect to qi instead are total derivatives in the sense that the dependence 
of p.(q), A(q) on q needs to be taken into account and the corresponding derivatives performed. The final 
result is then 


[Ds]i,j 


1 d'^n 


Vh? dqidqj 


A},/i,q=0 


(26) 
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The important point is that the above derivatives are taken for A;, /r constant even if they are themselves 
functions of q. Only the value of the order parameters at q = 0 is needed to calculate the superfluid weight. 
This is unsurprising since the superfluid weight is a linear response coefficient and therefore it depends 
only on the properties of the ground state. 

Supplementary Note 2: Wannier functions and tight-binding Hamiltonian 

Consider the Schrodinger equation with a periodic potential V{y) — V{y + a^), with the fundamental 
vectors that define the Bravais lattice, 

H'tp{Y) = + V (r)^ i/j(r) = £7p(r). (27) 

Bloch theorem states that the solutions can be labelled by a band index n and a quasimomentum index k 
and take the form^ 

V'nk(r) = , with eigenvalue e„k • (28) 

The functions gnk(^), called periodic Bloch functions (Bloch functions for brevity), are periodic in the real 
space coordinate gnk(j') = 5 nk(r + ai). The eigenvalues are periodic in the quasimomentum e„k = ^^(k+di)) 
with di the fundamental vectors of the reciprocal lattice that are defined by • dj = 27rSij. We adopt the 
convention that the periodic Bloch functions are normalized d^r |9rak(r)P = 1, where the symbol d^r 
denotes the integration over one unit cell. The band index n runs from 1 to +00 since there are an inhnite 
number of bands in the continuum. Consider now a subset S of relevant bands, also called composite band, 
separated from other bands below and above by sufficiently large bands gaps. It is possible to construct 
a particular convenient basis of localized functions that span exactly the same subspace spanned by the 
Bloch plane waves ipnk relative to the composite band. These are called Wannier functions and are defined 
by 

'tWa(l’- n) = f ^[C/k]a,nV'nk(l’) 

^ (29) 

'' ^ nes 

= ^ f 5][C/k]„,n5nk(r - n). 

Vh is the volume of the unit cell (in two dimension the area Aq of the unit cell is used instead) and the 
vector r; = ixSii + iya 2 + labelled by a triplet of integers i = {ix,iy,iz)'^, is a generic vector of the 
lattice. In the last line of the above equation we have made use of the fact that 5rtk(r) = 5nk(r ~ ri) 
for arbitrary i. The matrix I7k is a k-dependent unitary matrix that has to be chosen in order to obtain 
properly localized Wannier functions. 

The Wannier functions Wa{Y — r^) in Eq. (29) are obtained by the Fourier transform of the Bloch 
functions '*/’nk(r). The Fourier transform is a unitary operator, therefore the orthonormality of the Bloch 
functions is inherited by the Wannier functions 

y'^^i’V’*k(i’)V’n'k'(r) = ^^^<5n,nA(k-k') ^ J u;* (r - n)(r - Fj) = j . (30) 

According to Ref. 4 if the Chern numbers (one number in 2D and three numbers in 3D) of the composite 
band are zero then it is possible to choose C/k is such a way that the Wannier functions are exponentially 
localized. Moreover there are practical methods to calculate the matrices I7k that guarantee that the 
Wannier functions are maximally localized, more precisely they minimize a suitable localization functional 
(see Supplementary Note 5). The Wannier functions are useful since they are a basis of states for the 
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subspace of the Hilbert space relative to the given composite band and are very simple since they are 
translated copies by displacements ri of the small set of functions {?CQ-(r)}. The centers of the Wannier 
functions are defined as 

riQ, = + fya 2 + izUa + b„ , with ha = J (fr\wa{r)\'^r . (31) 

We now consider the Hamiltonian in second quantized form in the basis of Wannier functions. The braket 
notation is used for the Wannier functions (r|iQ;) = ^^(r —r;). The second quantized quadratic Hamiltonian 
is a — J (fir with ^jJ{r) the fermionic field operator satisfying the usual anticommutation 

relations {'i/'(r),= (5(r — r'). We expand the field operator in the basis of Wannier functions 

^(r) = + XI XI • ( 32 ) 

i aeS i 7 GS 

The functions w^{r — ri) are Wannier functions of the bands not included in the chosen composite band, 
and belong to the complementary set S. Inserting this expansion into the expression for Ti one obtains 


H — f d^r {r)H^{r) = X X4a(i«l^lj/3>%+ X 

a,d<SS ij 7,565 ij 

-EE Cia^iaJ/3Cj/3 + EE '^7-^i7J54i5 • 

a,/365 ij 7,565 ij 


(33) 


The key point in the above equation is that matrix elements of the form (ial H Ijy) between a Wannier 
function belonging to the composite band {a G S) and a Wannier function in the complementary set 
(7 G 5) vanish, the reason being that they are built from eigenfunctions of the Hamiltonian that belong 
to orthogonal subspaces. On the other hand off-diagonal matrix elements (ial H |j/3) with a, j3 G S are 
in general nonzero due to the fact that the Wannier functions are not eigenstates of the Hamiltonian. 
In particular there is a coupling between Wannier functions with a ^ jd due to the mixing given by the 
unitary matrix Hk in the definition (29). From the definition (33) it is also evident that the hopping matrix 
elements are translationally invariant ^fiaj /3 = Ka^jsii ~ j)- For exponentially localized Wannier function 
the matrix elements Ka^jsii — j) decay rapidly with the distance i — j between unit cells and it is usually a 
good approximation to retain only the matrix elements between a limited number of neighbouring Wannier 
functions. The result is the so-called tight-binding Hamiltonian^ which provides the physical picture of 
particles moving by discrete tunneling events between the sites of the lattice. 

Finally Eqs. (28)-(29) in the main text can be obtained by expanding the held operators in plane waves 
haa = Xlk {Nq is the number of unit cells) 


X X j/3Cj/3a = ^ J] X ( X 


a,^eS iJ 


a,peS i,j ^ ki 




k 2 


a,l3eS ki,k2 


(34) 


1 -J 


= E E4,7kdE'"'‘'*"‘''"'A'£.j« = E E4,[A'"(k)i.,«e„k, 

a,/365 k ^ i-j 2 k 


Supplementary Note 3: Superfluid weight in a multiband system at flnite temperature 

From now on we consider a translational invariant Hamiltonian and we use the notation employed in the 
main text. The symmetry constrain in Eq. (19) implies that for q = 0 the matrices £'k(q) and Wk(q) have 



the form shown in Eqs. (8) and (9) in the main text. It is convenient to rewrite the grand potential in 
Eq. (12) in a slightly different form which can be derived in a way similar to Eq. (22) 

fI(T,/., A,q) = ^ i(Tri7^(q) +TVi?^(-q)) - ^ ^ Tb [ln2(l + cosh/3Hk(q))] . (35) 

k k 


According to the results in the previous section the last term Tr|A(q)p/17 = can be 

dropped since it does not contribute to the superfluid density. Also the first term ^ (Tr i7^(q)+Tr i?^(—q)) = 
5 Xlk (Tr [£k-q ~ m1 ] + Tr [sk+q “ m1 ] ) does not contribute in the thermodynamic limit since is a pe¬ 
riodic function on the Brillouin zone. However we keep this term for reasons that will be clear in the 
following. Note that in the low temperature limit /3 —)• -|-oo the second term reduces to — jTr [|i?k(q)|] in 
agreement with Eq. (10) in the main text. 

The first derivative of the grand potential with respect to q is the current density 


J(q) = 


1 d^l 
Vh dq 


/i,A 


(Trl-^^k+q] - Tr[akek-q]) - ^Tr 

k ^ 


tanh 


j^ /?igk(q) ^ 


dqEi^{q) 


(36) 


The matrix elements of the diagonal matrix 9q£’k(q) can be calculated using the Hellman-Feynman theorem 
and the expression for i/k(q) in Eq- (6) in the main text 


[5qAk(q)]„_, 


W^(q)9qi7k(q)Wk(q) , 


with 5q^^k(q) 


^£k—q ^qT^k(q)\ 

-5qT>k(-q) ^kSk+q / ' 


We have defined X^k(q) = “^k-q^^k-i-q = ^k(“^) ™ the above equation. Notice that the Hellman- 
Feynman theorem holds only for the diagonal components a — b oi the matrices on the left and right hand 
side as emphasised by the square bracket notation [■]ab=a- The matrix tanh(£'k(q)) is itself diagonal and 
only the diagonal matrix elements in Eq. (37) are needed in the evaluation of the trace in Eq. (36). It is 
easy to check that the current is an odd function using the identities in Eq. (19). Again one can distinguish 
a conventional current contribution and a contribution that comes from the off-diagonal terms in 9qiFk(q) 
as in the zero temperature case discussed in the main text. Indeed the two contributions are separately 
gauge invariant as explained below. 

When taking an additional derivative with respect to q there will be terms proportional to the sec¬ 
ond derivative of the BdG Hamiltonian dq^dq.H]f^{q) and first derivatives of the unitary matrix 9qWk(q). 
The following identity is useful in order to express the off-diagonal matrix elements of Wk(q)9qWk(q) = 
— (<9qW^(q)) Wk(q) only in terms of the ground state solution 


Vy^(q)9qWk(q) 


Wi((q) aqiFk(q) Vyk(q) 


J a.h 


for a ^ b. 


^ cifi [Ak(q)]a,a [E'k(q)]b ,6 

The second derivatives of the quasiparticle excitaton energies can therefore be written as 

[dqMq^Eu{q)]a,a = [w^(q) 5g,ag.iFk(q) Wk(q)l 


(38) 


E 

b 

a^b 


1 


[E'k(q)]a,a - [Ak(q)] 6,6 


>T^(q) ^q^E\^{q) Vyk(q)J ^ ^ [W^(q) 5g,-ffk(q) Vyk(q)J ^ ^ ^ 


ih 


(39) 


Taking one more derivative of the current density and setting q = 0 leads to the superfluid weight 




\r 


^[dkAM - ^Tr 


tanh 


m 


dqA^Ek 


- -Tr 
2 


/3 


2 cosh"^(/3iJk/2) 


dq^Ekdq^Ek 


]} 


(40) 
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where i?k = £’k(q = 0), dq^E]^ = (9g.£’k(q = 0) and dq^dq^Ex^^ = ^q^^q^E\^{ci = 0). Eq. (40) can be written 
in a more convenient form by introducing some definitions (see Eqs. (8)-(9) in the main text) 


iVk,* = W^(q = 0) 5,.Hk(q = 0) Wk(q = 0) 

with 


_ _ [ -^k,i ^k,i 

y -^k,2 ^k,2y 

Ak,* = uldk^exJUx, + V^4,ekVk +Z^^9,,25kVk - V^^.^Pk^k 
5k,. = uldqPxMx. + vldq^Vx^Vx, + vldk^exMx. - uldk^ex^Vx, 


(41) 


[7k]a,6 = { 


/? 


2 cosh^(/3Ek/2)_ 

[tanh(/3Ek/2)]a,a - [tanh(/3Ek/2)]b^b 
[5k] a,a - [5k] 6,6 


for a = 


for a^b. 


(42) 


The abbreviation dq.'D]^ — 9g;X>k(q = 0) has just been employed and dq^dq^T>\^ = dq^dq.T>\^[f\ = 0) will be 
employed in the following. The hnal result for the superfiuid weight tensor contains three contributions 


[5s]i,i — 


1 


Vb? dqidqj 


— [5s,i]i,j T [5s,2]i,j + [5s,3]i,j j with 


[5: 




Z \—> 


/x,A,q=0 

Vk 


Vk+^^k ai^> dkA,£]i 


-PE> . . "k 


+1 


+ 1 


[D: 


s,2ji,i 


= — VT4 
vb? ^ 


Ux,Vl-Ux,- 


fiE> 

1 1 


+ 1 


-vi-Vk- 


e^^k + 1 




[5s,3]i,j — 2VK^ [[5k]a,6[5^k,.]a,6[5k,j]6,a ■ 


(43) 

(44) 

(45) 

(46) 


k a,b 


We have introduced the diagonal matrix E^ = diag(E„k), he. the upper diagonal block of Ek which 
contains the positive eigenvalues E„k > 0 (see Eq. (8) in the main text). The k-dependent effective mass 
tensor —3^— = ^ ^ enters the expression of Ds i which is the combination of the first term in 

^eff.k j ’ 

o 2 

Eq. (40) and the terms proportional to g ^ in the second term in the same equation (thus our derivation 

^ J 

holds also in the case of a finite system for which X]k'^[^k+q] is not a constant as a function of q). The 
analogous quantity that enters Ils ,2 is 

dq.dq^Vy, = {dk,Ql)i!S.{dkjQk) + idkjGi)A9kiSk) - (dkAj^k) ^^k - ^k^ {dkAj^^) ■ 

The above expression simplifies if the diagonal matrix A is proportional to the identity, resulting in (A is 
now a scalar, not a matrix) 


dqAAk = 2 A 


dkiGtdk,Gk + ^fci^k^^k 


(48) 


In deriving the above equation we have used the identity dkAjiGl^Qk) = 0, a consequence of the fact 
that Gx^Gk ~ 1. The component of the superfluid weight ^ results from the combination of the second 
term in Eq. (39) and the last term in Eq. (40). It contains both a term which is possibly hnite at zero 
temperature (see Eq. (14) in the main text) and a negative term which in the single-band case corresponds 
to the normal fluid component. I?s ,3 must be present in order to preserve the invariance of the superfluid 
weight, an observable quantity, under gauge transformations. In this context the gauge transformations 
are related to an inherent ambiguity in the definition of the Bloch wavefunctions namely the possibility 
of applying a transformation of the form Qx^ —> ^k^^k with ^k a unitary matrix subject to the additional 
constrain 

[ek,A]=0. (49) 
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This ensures that the form of the BdG Hamiltonian in Eq. (6) in the main text is preserved under 

the transformation. If the band dispersions are nondegenerate then is simply a diagonal matrix of 
phases = diag(e*'^"*^*^)). T>s,i is invariant under such a transformation since one has Vk —> »4,j^Vk and 
Eq. (49) implies that ^k commutes with the derivatives of e^. This means that also the sum Z)s ,2 + £*s ,3 
is gauge invariant. Below we show in a specific case that, by a suitable gauge transformation, it is possible 
to make the component Ds ,3 vanish leaving only Dg i and I?s, 2 - But it is not clear if this is possible 
in general and it is convenient to allow some freedom in the choice of gauge. The calculation of the 
superfluid weight in the multiband case requires only the solution of the ground state problem, namely the 
quantities Wk, Vk, Enk, A^, and the only additional ingredients with respect to the single-band case are the 
Bloch functions and their derivatives dk^Qk and dkidk^Qk that enter in Dgfi and I?s, 3 - It is easy to check 
that in the zero temperature limit Eqs. (41)-(46) reduce to the zero temperature result in the main text 
(Eqs. (12)-(14)). Eqs. (41)-(46) are very general and can be applied to a number of lattice models. 

We also provide explicit expressions in terms of the matrices Uk and Vk and the quasiparticle energies 
Enk of some correlation functions that are useful for the implementation of the self-consistency loop nu¬ 
merically or to check the self-consistency of an analytical solution. The first is the so-called anomalous 
Green function 


EiajjS — (Cia-l-Cj^l) 


_L p*k (riQ-rj,3) 

k 


Gk 



Uk- 


-VI 


Vk 


1 

+ 1 




(50) 


The unitary of Wk(q) implies that UkVl^ = VkUk ''^hich means that Eq. (50) is an Hermitian matrix 
Fiajp = E^p and, therefore, the pairing potential Aq, = U{ciai;Cia\) = UEia^ia is real. Another correlation 
function of interest is the usual Green function G 


G'lajp — (cf 


~ (4aiG/34) 


— V, 

Nc ^ 


-jk-(ria-rj;3) 


^k U^k 




+ 1 




+ Vk- 


,-PE> 


+ 1 


-V 


Sk 


a,j3 


(51) 


The total particle density is immediately obtained from the Green’s function riia = (cj^iCiQ-i-) -|- (c;« 4 ,G«t) = 

2GiQ.^iQ,. 


Supplementary Note 4: Superfluid density of the time-reversal invariant Harper- 
Hubbard model 

We now specialize our result for the superfluid density in a generic multiband system to the case of the 
Harper model. The band dispersions e„k and periodic Bloch functions gnk of the Harper model with 
commensurate flux per plaquette = 1/Q = 1 /Aorb are given by the Harper equation 


^ fi'nk(a) = -^5nk(a) , (52) 

with periodic boundary conditions gnk{c() = 5rak(<a-l-Q). A useful property of the Bloch functions that can 
be derived from the definition is 


E ^‘"gnkia -Fl) -h e * ^°-g„kia - 1) -F 2 cos ( k^a - 27r — 


a 


£ 

n 



Enk, f^n^k I k,^^ ~ 3nk{c>t 1) 


with 



(53) 


By means of the above identity it is shown below that a constant gap function Aq, = A is a self consistent 
solution of the gap equation, therefore in the following A is a positive real number and not a matrix as 
defined previously. This implies that Gj^^Gk ~ A^^^k = Al. In this case the BdG Hamiltonian i4k(q = 0) 
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can be easily diagonalized and the solution reads 


[^k]n,n^ ^nk.^n,n' 

[h^k;]n,n' “ Vn]i5n,n' 


with 

with 




'^nk — 


^/2 



En-k ) 


1 

2 


^nk M\^ 
En). ) ’ 


Enk = V {^nk - + A2 . 


(54) 

(55) 

(56) 


A fundamental object for what follows is the projector operator P„k defined in terms of its matrix elements 

[Pnk\a,P — [Qk\^^JlQ]^n^p- (57) 

This is a projection operator since it is positive P„k > 0 and idempotent = P„k- The above definition 
is immediately generalized to the projection operator onto a given set of bands simply by summing over 
the corresponding subset of values of n. Note that Eq. (53) imply that 






= 1 . 


(58) 


Thus the gap equation at zero temperature reads 


^ — UE\gi i(Ti — ^ ^ Unk^nk tanh [-EnkJojC 


k.n 


wX i: 


7 ^ , , PEnk 

, , -tanh- > 

Nc^ ^ 2Enk ^ 

n kSred.B.Z. 


P ^ ^ A 

^ 2Enk 

^ n kSred.B.Z. 


tanh 


2 

7 

PEnk 

2 ’ 


(k+^k.) 


(59) 


In the second and third line of the above equation the vector k takes values on the reduced Brillouin zone 
kx,ky e [—7r/(Qa),vr/(Qa)] since the degeneracy due to Eq. (53) has been taken care of by the summation 
over m in the second line. Analogously the average number of particles per lattice site rip = riia is fixed 
by the equation 


Up 

2 


E E \ 

n kGred.B.Z 



^nk M , 1 

—— tanh 


PEnk \^ 


(60) 


If a solution of the coupled equations (59)-(60) is found then our initial choice of a constant pairing potential 
is self-consistent and this is in general the case since we have two equations with two unknowns. 

Up to now all of our calculations have not involved any approximations other than the mean-field 
one. We now use some approximate results for the eigenfunctions and eigenvalues of the Harper model 
as detailed in Ref. 5. It has been shown that for Q ^ I the first lowest laying bands are separated by 
a gap which is an algebraic function of Q, while the bandwidth of each band is exponentially suppressed 
(max^enk ~ miiikEnk oc 6““*^ with awl). This corresponds to the the fact that the Harper model bands 
approach the perfectly flat Landau levels in the continuum limit. We call £n ~ Enk the average energy 
of the h-th band and choose a value of the interaction strength U such that Eq. (15) in the main text is 
satisfied, meaning that interband effects and effects due to a nonflat band dispersion are neglected. As a 
consequence the summation over the wavevector k produces N^/Q = NcU^ identical terms in Eq. (59)-(60). 
We write np/{2n^) = h -|- 1 ® with h the integer number of the highest (partially or completely) occupied 
band of the Harper model and with 0 < v < 1 the filling of factor of the same band. The filling factor i/ is 
spin resolved and is the same for up and down spins. The bands are numbered starting from n = 0 up to 
n = Q — 1 to conform with the usual Landau level notation. In the flat band approximation e^jk ~ £«) one 
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has from Eq. (56) that Efik ~ En and since we expect U ^ En En]i for n n, only the n = n term is 
relevant in Eq. (59). By cancelling A on both sides one obtains the equation for En 


a = Ih* tanh ^ 


(61) 


This equation has nonzero solution only for /3Un^/4 > 1. Indeed ksTc = Un^/A is the mean-field critical 
temperature at half filling (see below). Close to the critical temperature it is possible to use the approxi¬ 
mation tanhx ss X — ^ and derive En{T < Tc) w 2\/3kBTy/l — T/T^. Using Eq. (60) as well it is easy to 
obtain the approximate self consistent solution at finite temperature 


En, n = n 

Jenk-^l, nj^n, 

1 


Enk 


fi = en + Un^[iy - - 


A = 


Un^ / ( 

2 \l\Un^J 




(62) 

(63) 

(64) 


If the argument in the square root of Eq. (64) is negative then A = 0 even if Eq. (61) has a solution. 
This can happen away from half-filling. Instead at half filling A = En. The zero temperature limit of 
the above solution is given by Eq. (16)-(19) in the main text. Since the Harper model has been recently 
implemented in ultracold gases we test the validity of the isolated flat-band approximation (Eqs. (61)-(64)) 
against the full self-consistent solution of Eq. (59)-(60) in Supplementary Figure 2. For reasonable (not 
too small values) of and U the flat-band limit is indeed realized to a good approximation. 

In the derivation of the approximate self-consistent solution (61)-(64) we have used the property (53) 
specific of the Harper model. However it is generally valid in the case where Aq, = A with A a real scalar. 
Indeed simply using this assumptions, which leads to Eqs. (54)-(56), and the condition (15) in the main 
text one obtains the two coupled self-consistency equations 



The projector on a single band is normalized according to TrFjjk = 1 and this implies N~^ '^k.[Enk]a,a = 
N~^. The result in Eqs. (61)-(64) is thus recovered with = N~^. 

We are now in position to evaluate the various contributions to the superfluid weight. The diago¬ 
nal superfluid weight is zero in the flat-band approximation. The other contributions can be evaluated 
straightforwardly and the result is 

2 A^m 

(67) 

where the temperature enters only in the energy gap A(r). This result reduces to Eq. (20) in the main 
text at zero temperature. In deriving the above equation we have used the fact that Unk^nk ~ 0 for n ^ n. 
We have also neglected terms with n h in Eq. (46) since En En^n according to Eq. (62). The matrix 
= Re(A4ij) is the real part of a Hermitian matrix defined as 


Mij 


27r 

'a 


E%(k) 


1 

27r 





1 

271 






+ Tr 


0l(dkA)^l(dk^Qk) 


( 68 ) 
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The integration is over the full Brillouin zone. The hrst term in parenthesis in the above equation corre¬ 
sponds to Ds ,2 while the other one corresponds to Ds, 3 - Crucially the matrix A4 is defined in terms of the 
matrix which is now a Q x 1 column vector (in general a rectangular matrix), defined as the projection 
of Qk on the Bloch wavefunctions corresponding to the h-th band 

[Q\i\a,l = [Qk]a,n , = 1 , GkGk ^ • (69) 


It can be easily shown that if Qk is a square unitary then Ai is zero. In our context in fact a square unitary 
matrix is a pure gauge transformation which must not have observable consequences. The matrix A4 is 
also positive semidefinite > 0 since it can be written in the form 


A4ij 


1 

27r 


d^k2Tr 


'B.Z. 


{dkiGk){^ - Pnk)idkjGk) 


(70) 


The invariant quantity I3ij(k) is called the quantum geometric tensor, its real part is the quantum metric, 
while the imaginary part of Bij(k.) is the well-known Berry curvature, and its integral over the Brillouin 
zone is the Chern number Im(AIij) = A4]j = cijC, generalized to a set of bands (a composite band)^ 
(cij = —Cji is the Levi-Civita symbol). The second term in parenthesis in Eq. (68) is real and does not 
show up in the Berry curvature. In can be easily verified that the quantum geometric tensor is an invariant 
under gauge transformations Gk Gk^k where is a unitary matrix with dimension equal to the number 
of columns of Gk- 

We now want to evaluate Aly using a suitable approximation for the wavefunctions of the Harper model. 
Such an approximation can be obtained by taking the Harper equation (52) and replacing the cosine 
potential with its second order expansion around the minimum and furthermore replacing the discrete 
differences with derivatives. We obtain the equation of an ordinary harmonic oscillator for the Bloch 
wavefunction V’(riQ,) = Qiy) 

- V'(a) = with ^ • (71) 


The variable a is now a continuous variable. We are neglecting effects due to the tunneling between the 
minima of the cosine potential in the Harper equation (52) in the main text. These corrections can be 
shown to be of the same order of the exponentially small bandwidth® which we neglect in our treatment. 
The above equation is precisely the equation that one obtains when solving the problem of a quantum 
mechanical particle in a magnetic field introduced by means of the Landau gauge (note the shift of the 
harmonic potential proportional to kx)- The solutions are the harmonic oscillator wavefunctions iPn{x) 
(Gaussian function times an Hermite polynomial) which corresponds to the different Landau levels. It 
is also useful to introduce annihilation and creation operators to manipulate such functions (x — x and 

p = -idx) 


aipn = \/npn-i , (72) 

aVn = Vn + l(pn+i ■ (73) 

Eq. (71) can be rewritten in terms of the annihilation and creation operators in the form 

^)V'(a) = ^(j+4)V'(a)- (74) 

The lattice periodicity is taken into account by imposing that the Bloch wavefunction transform in the 
proper way under translations, namely one has 

V-„k(ria) = E (a + Qiy - Qs - . (75) 
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The above wavefunctions satisfy the properties required by Bloch theorem, namely V'nk(ria + ax) = 
'i/'nk(ria + aQf) = e*'^^!'“V'nk(ria)- In Order to ease the notation in the following we 
label the positions of the lattice sites in the square lattice in the more natural way r; = a{ix^ + iyy)t rather 
than Tic, = a{ix^ + (a + Qiy)y)- Checking that the above wavefunctions are orthogonal and normalized is 
a good way to get a better feeling of the approximations that have been made (the number of lattice sites 
is Ng = L? with L = QM and the therefore then number of unit cells is Nc = Ng/Q = QM^) 


J{kx-K)jxa iQ{kyS-k'yS')a, 


I]V’„k'(rj)V'nk(rj) = 2 

1 Q 


S,S Jx,Jy 


X ‘Pn (^jy -Qs- Pn (^jy “ Qs' - 

- E E (jy -Qs- Pn (Jy - Qs' - 


M 

^kxK 

M 


S^S' jy 


27r 


(76) 


Y e^Qiky-k'y)sa^2^ E - Qs - 9^\ = ^ ^iQ{ky-k'y)sa 

S jy ^ '' S 


= hx,k'jky,k'y ■ 


We have neglected terms where s ^ s' since the harmonic oscillator wavefunctions are well localized inside 
the large {Q sites) magnetic unit cell and their overlap with wavefunctions in neighbouring magnetic unit 
cell is exponentially small. Moreover we approximate the sum over jy with the integral ’YhjyP'nijy) ~ 
j dx Pnix) = 1, which leads to an error of the same order as it can be checked with the Euler-MacLaurin 
formula. The harmonic oscillator wavefunctions are normalized in the usual way. The periodic Bloch 
functions gnk{jy) are thus (the same as Eq. (25) in the main text with a different labelling) 


9nUjy) = E Vn (i, - Qs - . 


(77) 


We can now calculate Sij(k) within this approximation for the periodic Bloch functions. We consider hrst 
the quantity Gl.{dkyGk) that enters the second term in parenthesis in Eq. (68) 

Q 

^k(^fca:^k) = 9nk{jy)^kj:9nk{jy) 
jy — ^ 

= E (i, - <3» - 0 t. 9 n (i, - <3*' - 

= E E »>»(*- (* - e* - ^) 

jy = l S ^ 7 \ / 

" E (7;/) J^{a - d'f)ipnijy) 

jy 

= JI E (jy) (vn(pn-i{jy) - Vn + lpn+i{jyfj = 0 . 

3y 

We have combined the sum over jy = 1,... ,Q and the sum over s in a single sum over all integer jy. 
This sum over jy has then been approximated by an integral and the overlap of orthogonal functions 
ipn carried out as usual leading to a zero result. One can check that the choice of the periodic Bloch 
functions in Eq. (77) leads to Gl^(dkyGk) 7^ 0- However we can perform a gauge transformation of the 
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form gnk —)• Going through the calculation along the same line of Eq. (78) one obtains 


something proportional to Ylijy 3y^‘n{jy) ~ / dxxip%{x) = 0. We will make use of this gauge choice when 
convenient. 


First we calculate Bxyi}^)- We only need to evaluate Tr \^dk^Q]^{dkyQ\!^ 
parenthesis in Eq. ( 68 ) is zero according to Eq. (78) 


since the second term in 


Q 

Bxyi}i-) ~ Tr (dk^Q\^{dkyQ'k) ~ ^k^9nk{jy)^ky9n]i{jy) 

J!/ = l 
Q 


jy = l S,S' 


e^ky{3y-Qs)a^_ _ 


QhxCi 

2tt 


dky 


Y^dk 


Jy 

iQa? 

iQa^ 

iQa? 

All 


^ikyjya^^_ 


'Pn 1 Jy 


Qkxd 

2tt 


dk 


Pn 1 3y 


^_,ky[jy-Qs')a^_ _ g^/ _ 


Qkxd 

27r 


/ Q kxd 

27r 


E , (. Qkxa\ . (. Qkxd\ 

Pn[3y 2 ^ j3yPn[3y j 

]y 


(79) 


3y 


T. 


1 Jy 


Qkxd 

2tt 




Qkxd 

2-k 


Pn 1 3y 


Y. ^Pn-l{3y) -n + lifn+lQy) Vnpn-lijy) + Vn + lipn+l{jy) 


Qkxd 

27r 

iQd? 

dvr 


3y 


Next we evaluate Bxx()^)- Again we need only need to evaluate Tr \ {dk:^Q]Q{dk^Q\^)^ according to Eq. (78) 
1 ^ « 


^Sxx(k) — Tr ^k:r39n]<.{jy)^kx9n\<i{3y) 

jy — ^ 

E o / . QkxCL\ f . Qkx(l\ 

dk^Pn J dk^Pn J 

) Pn (jy) 


Qa? 


E 

Jy 


LO 


UJ 

2 


) Pn {jy) 


(80) 


dvr 


= Y 'dnpn-l{jy) - Vn + lipn+l{jy) Vnipn-l{jy) “ VU + lipn+l{jy) 


Jy 


(2n+l). 


For Byy(k) we choose the modified gauge that ensures the vanishing of the second term in parenthesis in 


16 



Eq. (68) 


^Syy(k) — Tr 9kygnk.{jy)dkygn\<i{jy) 

3y — ^ 


= E«‘. 

3y 


3y 


Jy 


■ _QkmR\n 
y ) ifin 

(^jy- 

Qkxd 

2n 

Qkxd^ 1 

3y- 

Q kxa 

2n r^\ 

27r 




e Miy ,jy 


QkxCi 


271 


Jy 


Qkxa 

27T 


\ Jy 


Qkxa 

27T 




Qa? 

An 

Qa? 

An 


\f^ 


Vnipn-i{jy) + Vn + lipn+i{3y) Vn^Pn-iijy) + Vn+ l(pn+i{jy) 


3y 


(2 n+l). 


(81) 


The quantum geometric tensor Bij is constant over the Brillouin zone and the calculation of A4 is trivial 
(the integration is over the full Brillouin zone kx — [—tt/o,tt/o], kx = [—n/{Qa),n/{Qa)\) 

Qa^ 


>, = i/d^kB(k) = ^ 


dk 


/B.Z. 


2n + 1 —i 
i 2h + 1 


2n + 1 —i 
i 2h + 1 


(82) 


We obtain the result that for a Landau level the Chern number is |(7| = 1 and that the superfluid density 
in the h-th Landau level is proportional to 2h + 1. In particular for the lowest Landau level the bound in 
Eq. (23) in the main text is saturated. 


Supplementary Note 5: Bound on the localization functional for Wannier functions 


Marzari and Vanderbilt consider the following functional that can be minimized numerically in the matrix 
C/k ill Eq. (29) in order to construct maximally localized Wannier functions® (we consider the 2D case for 
simplicity) 


E = |^(Oa| lOa) — |(Oq;| r |Oa)p , 

a 


(83) 


(Oa| |Oa) = / d^r |u;a(r)|V^ =/ d^k (Vk5ak|Vk5ak) , (84) 

Jn (271-)^ Jb.z. 

(iQ;|r|0^)=y (frWa{r-ri)wi3{r)r = ij0^J d^ke*’^'*''(5rak|Vk5/3k) ■ (85) 

Aq is the area of the unit cell, and the periodic Bloch functions indexed by a instead of n are obtained by 
the transformation = I]„[t/k]a,n5nk- We have also dehned the ket |Vk5ak) as (r|Vk5ak) = Vk5ak(r) 
in the identities (84)-(85) which are derived in Ref. 6. The above functional can obviously be interpreted 
as the sum of the spreads of the Wannier functions. The localization functional can be split into two 
contributions 



F = Fi + F, 

(86) 

7---E 

(Oa Oa) — (i/? r Oa) p 

(87) 

a 




F= Y1 Ki/3|r|Oa)y 

(88) 


id 

i/S^Oa 
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The first contribution Fj is gauge invariant in the sense that it is independent of the unitary matrix 
C/k and it is interesting to provide an expression for Fj in terms of the Bloch functions Qak using the 
identities (84)-(85) 


i a 

-E 


An 

( 2^)2 

( 2 vr)' 


B.Z. 

2 




'B.Z. 


(fkl[ (4i5ak2l5/3k2> (5/3ki|4.5aki) 

Jb.z. 


^ ^ ^ ^ ^ y {9ki9ak\dki9ak) + ^ ^ (5«k|^fci5/3k) {9fik\dki9ak) 

\ ''J ,• JB.Z. n ft 


(89) 


An 

(27r)2 

An 

w? 


Y, / [Th[(aial)(a,gk)] + {d^Gu) Gl (ai^k)]] 


I B.Z. 


Y[ d\lBu{k)^ 

, Jb.z. ^ 


1 

^ 2 


Tr[A4] 


In the last line we have used the definition of quantum geometric tensor Bij (k) and of the invariant matrix 
A4 [see Eq. (68)]. In particular the matrix [Gk]r,a ~ 5 ak(r) is analogous to the matrix Gk as defined 
in Eq. (69) and used in the main text with the only difference that the index r is continuous and has 
to been integrated on the whole unit cell when the matrix product is performed. Using the arithmetic 
mean-geometric mean inequality (Ai -f A2)/2 > \/AiA 2 for the eigenvalues of the positive matrix A4 and 
indicating with the real part of M. and with the imaginary part, one has the following chain of 
equalities / inequalities 

iTr[A4] = ^Tr[A4^] > ^det[A4R] > ^Jdet[JkF] = \C\. (90) 

The second inequality is Eq. (23) in the main text. Thus the Chern number of the composite band provides 
a bound from below of the localization functional, namely the result stated in the main text 

E>^|C|. (91) 


Supplementary Note 6: Berezinsky-Kosterlitz-Thouless transition 

The low energy fluctuations of the order parameter phase 6{y) have an energy cost given by the Hamiltonian 
of the classical 2D XY model 

nxY = ^ J d^v{V9{r))\ (92) 

The “stiffness” T is in fact related to the superfluid weight by the relation T = Dg/i^/d where the factor 4 
that appears in this last equation is due to the fact that in our notation V0(r) = 2q for a constant phase 
gradient. At the Berezinsky-Kosterlitz-Thouless transition temperature Tbkt the vortex unbinding due 
to thermal fluctuations leads to the destruction of the phase order and thus of superfluidity. The BKT 
critical temperature can be obtained by the universal relation^ E'(rBKT)/(fcB7BKT) = 2/7r. We would like 
to estimate the BKT critical temperature with the hnite temperature result obtained for the superfluid 
weight (67). Eor simplicity only the half hlling case ^ = 1/2 is considered. The zero temperature value of 
the energy gap A{T = 0) = Un^/2 = 2kBTc is taken as the energy scale. The temperature-dependence of 
the energy gap is then obtained by the equation (compare with Eq. (61) and (64)) 


A(r) 

A(0) 


= tanh 


fTcA{T)\ 

\T A{0)) ’ 


(93) 
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while the stiffness is, from Eq. (67) 


TniT) _ 2h + l A2(T) 

A(0) ~ 47r A2(0) ■ 

The dependence of the stiffness on the Landau level number n has been introduced 
that needs to be solved in order to find the BKT critical temperature is 

2h + 1 A^(Tbkt) _ Tbkt 
^ A2(0) Tc ■ 

Eq. (93) and (95) can be solved numerically and the result is shown in Supplementary Figure 1. The BKT 
temperatures for h = 0, 1 , 2 are approximately Tbkt ~ 0.25, 0.61, 0.75 Tc, respectively. This justifies our 
assertion in the main text that the mean-field temperature is in fact a good estimate of the real temperature 
of the superconducting transition in the lowest bands of the Harper model. 


(94) 

. Then the equation 
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